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We argue that magnetic flux-conservation in turbulent plasmas at high magnetic Reynolds numbers nei- 
ther holds in the conventional sense nor is entirely broken, but instead is valid in a novel statistical sense 
associated to the "spontaneous stochasticity" of Lagrangian particle trajectories. The latter phenomenon 
is due to the explosive separation of particles undergoing turbulent Richardson diffusion, which leads to a 
breakdown of Laplacian determinism for classical dynamics. We discuss empirical evidence for spontaneous 
stochasticity, including our own new numerical results. We then use a Lagrangian path-integral approach to 
establish stochastic flux-freezing for resistive hydromagnetic equations and to argue, based on the properties 
of Richardson diffusion, that flux-conservation must remain stochastic at infinite magnetic Reynolds number. 
As an important application of these results we consider the kinematic, fluctuation dynamo in non-helical, 
incompressible turbulence at unit magnetic Prandtl number. We present results on the Lagrangian dynamo 
mechanisms by a stochastic particle method which demonstrate a strong similarity between the Pr m = 1 and 
Pr m — dynamos. Stochasticity of field-line motion is an essential ingredient of both. We finally consider 
briefly some consequences for nonlinear MHD turbulence, dynamo and reconnection. 
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I. INTRODUCTION 



Hannes Alfven in a seminal paper in 1942 introduced 
the notion of flux-freezing in magnetohydrodynamic plas- 
mas at infinite conductivity, noting that "every motion 
(perpendicular to the field) of the liquid in relation to 
the lines of force is forbidden because it would give in- 
finite eddy currents'^. In the years since, the property 
of flux- conservation has become a powerful tool in the 
analysis of many near-ideal plasma phenomena. For ex- 
ample, in his excellent monograpfP, Kulsrud states that 
"The most important property of an ideal plasma is flux- 
freezing" , before proceeding to illustrate its many appli- 
cations. Of course, physical plasmas in the the labora- 
tory and in astrophysics are subject to various forms of 
non-ideality, including Spitzer resistivity, ambipolar dif- 
fusion, etc. The general assumption in the field of plasma 
physics, however, is that as long as such non-ideality is 
sufficiently "small" , then flux-freezing will hold in an ap- 
proximate sense. This idea dominates the discussion of 
turbulent magnetic dynamo at high kinetic and magnetic 
Reynolds numbers. For example, it is commonplace to 
find statements in the literature such as the following: 
"The small-scale turbulent dynamo is caused by the ran- 
dom stretching of the (nearly) frozen-in field lines by the 
ambient random flow" 3 , or "It is well established both an- 
alytically and numerically that a weak magnetic field can 
be amplified by the random motions of a highly conduct- 
ing fluid [1-3]. This occurs because magnetic-field lines 
are generically stretched by the random motions of the 
fluid in which they are (almost) 'frozen.' "S. The phys- 
ical idea which underlies these statements is that a tiny 
resistivity should diffuse field-lines only a short distance 
through the plasma. For example, the usual quantitative 



estimate is well expressed in this quote from Kulsrud's 
monograph, Ch.13, on magnetic reconnection: 

"Flux freezing is a very strong constraint 
on the behavior of magnetic fields in astro- 
physics. As we show in chapter 3, this implies 
that lines do not break and their topology is 
preserved. The condition for flux freezing can 
be formulated as follows: In a time t, a line of 
force can slip through the plasma a distance 



(1) 



If this distance £ is small compared to 5, the 
scale of interest, then flux freezing holds to a 
good degree of approximation." 

We shall argue that these commonplace ideas on flux- 
freezing are wrong. They contain an implicit assumption 
that the plasma fluid remains smooth and laminar for 
very small non-ideality. The quantitative estimate that 
field-lines slip through a resistive plasma only a diffusive 
distance cx y/rjt in time t is incorrect — by many orders of 
magnitude — in a turbulent plasma. Since laminar flow at 
very high kinetic and magnetic Reynolds is unstable to 
development of turbulence, flux-conservation in the con- 
ventional sense must be the exception rather than the 
rule in astrophysical plasmas. Indeed, we shall show that 
the standard views on flux-freezing must be incorrect, be- 
cause the very notion of a Lagrangian fluid particle tra- 
jectory breaks down in turbulent flow with a spatially 
"rough" velocity field (i.e. with a power-law kinetic en- 
ergy spectrum similar to that of Kolmogorov). Recent 
research has discovered a novel phenomenon of "spon- 
taneous stochasticity", according to which fluid particle 
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trajectories are intrinsically random in high-Reynolds- 
number turbulenc e 5 " 12 !. This surprising phenomenon is a 
long overlooked consequence of the fluid-dynamical effect 
of Richardson two-particle turbulent dispersion^ Be- 
cause of spontaneous stochasticity, it makes no sense to 
assume that a field line follows "the" plasma fluid ele- 
ment, because there are infinitely many distinct fluid tra- 
jectories starting from the same point! But it is also not 
true that flux-freezing is completely broken. We shall ar- 
gue below that magnetic flux-conservation remains valid 
in the ideal limit of high Reynolds numbers, but in a novel 
stochastic sense associated with the intrinsic stochastic- 
ity of the Lagrangian particle trajectories. 

A correct formulation of flux-freezing is fundamental 
to understand a number of important astrophysical pro- 
cesses, such a s tur bulent dynamo and reconnection. In 
previous worlfHEsl we have shown how stochastic flux- 
freezing is involved in the small-scale "fluctuation dy- 
namo" for a soluble model problem: magnetic fields ad- 
vected by the Kazantsev-Kraichnan ensemble of veloc- 
ity fields that are spatially rough and white-noise in 
time^HUl. it is worth noting, by the way, that "spon- 
taneous stochasticity" is a rigorously established phe- 
nomenon for this modePa i t 

was shown that presence 
of fluctuation dynamo effect at zero magnetic Prandtl 
number depends crucially on the degree of angular cor- 
relation between the infinite number of magnetic field 
vectors that are simulta neous ly advected by turbulence 
to the same spatial pointP3E33. 

Here we shall make a similar study for kinematic 
dynamo in non-helical, hydrodynamic turbulence, by 
a Lagrangian numerical method that employs data 
from a hig h Reynolds-number turbulent flow archived 
online^Ml \y e present results for unit magnetic Prandtl 
number which demonstrate — and quantify — the effect of 
Richardson diffusion and stochastic flux-freezing on the 
small-scale turbulent dynamo. We find, in fact, remark- 
able similarities between the Lagrangian mechanisms of 
small-scale dynamo in hydrodynamic turbulence at unit 
Prandtl number and in the Kazantsev m odel at zero 
Prandtl number. Previous numerical studies^!! have in- 
stead found a close similarity of the unit Prandtl-number 
fluctuation dynamo with the solution of the Kazantsev 
model at infinite Prandtl-number, observing, in particu- 
lar, the "Kazantsev spectrum" fc 3 / 2 of magnetic energy 
at high wavenumbers. One studj0went so far as to claim 
that " It is clear that the kinematic dynamo in such 
[Pr m = 1] runs is of the large- Pr m kind." This inter- 
pretation is ruled out by our new results at much higher 
Reynolds numbers, which show that the inertial-range 
phenomenon of Richardson diffusion strongly affects the 
exponential growth rate of magnetic energy in the kine- 
matic regime. However, stochastic flux-freezing is not a 
property of kinematic dynamo only but will hold also for 
fully nonlinear MHD turbulence and have important im- 
plications there for magnetic dynamo and reconnection. 

The detailed contents of this paper are as follows: In 
the following section II we shall briefly review the phe- 



nomenon of "spontaneous stochasticity", both its theo- 
retical bases and its present confirmation from simula- 
tions and experiments. We aso present new numerical 
results of our own which support the essential predic- 
tions. In Section III we discuss stochastic flux-freezing. 
We begin with a demonstration of the stochastic flux- 
conservation properties of resistive MHD, which are novel 
results themselves. We employ Lagrangian path-integral 
methods that provide good physical insight. We then dis- 
cuss the ideal case, via zero-resistivity and other limits. 
In section IV we employ the new results to discuss the 
turbulent kinematic dynamo. After reviewing the La- 
grangian theory of dynamo, we present our numerical re- 
sults and their comparison with analytical results for the 
Kazantsev model at zero Prandtl number and with pre- 
vious numerical studies. In section V we briefly discuss 
some implications and open problems for nonlinear MHD 
turbulence and section VI contains our final discussion. 
An Appendix sketches the derivation of the path-integral 
formulas used in the main text. 



II. RICHARDSON DIFFUSION AND SPONTANEOUS 
STOCHASTICITY 

A. Richardson 2-Particle Dispersion 

We briefly review Richardson's theory^ of 2-particle or 
relative turbulent dispersion, emphasizing perspectives of 
recent research. See also more complete review a 22 ™ 24 !. 

The object of Richardson's study was the separation 
Ax(t) = xi(i) — X2(i) between a pair of passive La- 
grangian tracer particles in a turbulent flow, such as 
ash particles in a volcanic plume. Richardson's approach 
was semi-empirical. By estimating the "effective diffu- 
sivity" K = (|Ax| 2 )/i as a function of rms separation 
£ = yj (|Ax| 2 ), he inferred from data that there is a scale- 
dependent diffusivity coefficient 

K(£) ~ (1) 

This is essentially a "running coupling constant" in the 
modern sense of renormalization group theory. Richard- 
son proposed further that the probability density func- 
tion of the separation vector £ = xi — X2 would satisfy a 
diffusion equation 

d t P(£,t) = ^(K(£)^(£,t)y (2) 

Richardson observed that there is an exact similarity so- 
lution of his equation given by a stretched-exponential 
PDF, which he wrote explicitly in one space dimenson. 
Here we note its form 

in the more physically relevant case of three space di- 
mensions. All solutions of |2]) approach this self-similar 
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form asymptotically at long times^. Averaging £ 2 with 
respect to the self-similar density ^ yields 



(t 2 (t))= 7o t 3 



(4) 



with 70 = ^^Kq . This is the famous Richardson £ 3 -law. 

Richardson's work preceded the Kolmogorov 1941 
(K41) theory of turbulence, but it was shown by 
Obukhov^ to be fully consistent with that theory. This 
can be seen by a toy calculation in one space dimension. 
Assume that l[t) satisfies the initial-value problem 

with velocity increment 8u(£) scaling as in K41 theory, 
where e is the mean energy dissipation per unit mass. 
Separation of variables gives the exact solution 



*(*) = el /3 + (goe) 1/3 t 



3/2 



(5) 



2 /3 

if one defines a time t$ = £ Q 
izes the initial separation then, for t 3> to, 



/(goe) 1 ^ 3 which character- 



e 2 (t) 



g et 3 . 



(6) 



For sufficiently long times, the particles "forget" their ini- 
tial separation and the Richardson law is obtained. The 
dimensionless parameter go which appears in this form of 
the f 3 -law is usually called the Richardson- Obukhov con- 
stant. The physical mechanism of the explosive separa- 
tion of particles, even faster than ballistic, is the relative 
advection of the pairs by larger, more energetic eddies as 
their separation distance increases. 

This physics seems relatively simple and benign, but it 
has extraordinary consequences. As first pointed out in 
a seminal paper of Bernard, Gaw§dzki, and KupiainerP, 
Richardson's theory implies a breakdown in the usual 
notion of Laplacian determinism for classical dynamics! 
This may already be seen in our toy calculation above. If 
we set Iq = 0, then the solution ^ becomes the Richard- 
son law i 2 (t) = go£t 3 > for all positive times t. Thus, 
two particles started at the same point at time sep- 
arate to a finite distance at any time t > 0. The same 
oddity may be seen in Richardson's similarity solution 
([3]), which satisfies at initial time t = 

P,(£,0)=d 3 (£). 

All particles start with separation £(0) = 0. However, 
P*{t,t) is a smooth density for t > 0, so that £{t) > 
with probability one at later times. Richardson's the- 
ory thus implies that two particles advected by the fluid 
velocity u(x, t) which start at the same initial point Xq 



dt 



x(i)=u(x(t),i), x(0)=x 



can follow different trajectories. This seems to violate the 
theorem on uniqueness of solutions of initial- value prob- 
lems for ODE's. However, such theorems assume that the 



advecting velocity u(x, t) is Holder- Lipschitz continuous 
in the space variable x 



|u(xi,t)-u(x a ,i)| < C|xi-x 2 | 



(7) 



with exponent h > 1. A turbulent velocity field in a 
Kolmogorov inertial range has instead Holder exponent 
h = 1/3 and the uniqueness theorem need not apply. Our 
toy calculation earlier is just the standard textbook ex- 
ample for failure of uniqueness (see HartmarP^, p. 2) . In 
that example, for any non-negative "waiting time" r > 



^) = ( 50 e) 1/2 (t-r) 



3/2 



0. 



is a solution of the initial-value problem with £q 
(Here (x) + = x for x > and = otherwise.) 

The above considerations may seem fairly technical 
and mathematical. It has been shown, however, that 
this breakdown in uniqueness of trajectories can appear 
in various physical limits for turbulent advection. Even 
more remarkably, the solutions of the deterministic clas- 
sical dynamics become int rinsica lly stochastic! See the 
important series of papers^ 7 * 9 * 10 !. Advanced pro bablis tic 
techniques have obtained the most refined results^^. A 
very clear and concise review of the subject is available, 
which the reader may consult for further details. Below 
we briefly describe the most basic theory and results. 



B. High-Reynolds-Number Limit and Spontaneous 
Stochasticity 

The easiest way to understand the phenomenon is via 
the problem of stochastic particle advection, 



dt 



x(t) = u l/ (x(t),i) + \/2^r 7 (i), x(t )=x (8) 



with advecting velocity perturbed by a Gaussian white- 
noise rj(t) multiplied by V2k and with velocity assumed 
spatially smooth at subviscous length-scales I < l v , for 
a finite viscosity v. The transition probability for a single 
particle in a fixed (non-random) velocity realization \i v 
can be writte n using a "sum-over-histories" approach as 
a path-integra l 5 ' 28 1 29 1 : 



GS'"(x / ,t / |x 0j to)= / Dx^xj-xCt,)) 

Jx(i )=x 

x exp J* dr |x(t) - u"(x(r), r)| 2 ) . (9) 



Since this formula plays an important role in our analysis, 
we provide a self-contained derivation in the Appendix. 
A physical motivation to study such random advection 
is the problem of the evolution of a passive scalar, such 
as a temperature field or dye concentration. These fields 
solve the scalar advection-diffusion equation 



d t 0+(u u -V)0 = kA9, 



(10) 
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with k the molecular diffusivity. The ex act soluti on of 
( 10 1 is given by the Feynman-Kac formula^ 1 22 * 28 1 29 1 : 



x,i)= / d 3 a;o0(xo^o)G^ K (x o ,to|x^) 



a(t)=x 



Va6(a(t ),t 



xexp(-i-^ dr|a(T)-u>(T),r)| 2 ) (11) 

for to < t. This corresponds to solving backward in time 
the stochastic equation 



a(r) = u 1/ (a(r),r) + V2K^(r) 



from r = t to t = to, with the condition &(t) = x. The 
present value of the scalar field is thus the average, along 
stochastic Lagrangian paths, of its earlier values. 

It naively appears by an application of the Laplace 
asymptotic method to ^ that the transition probability 
collapses to a delta-function 

Gr(x/ I t/|x ,to)^5 3 (x / -x(t / )) (12) 

as k — y 0, with x(t) the solution of the ODE x = u(x, t) 
for initial condition x(to) = xo. Only for such time- 
histories is the action vanishing in the exponent of the 
path-integral. However, it was showiPH that (pi may 
not hold if simultaneously v — > (or the Reynolds- 
number Re = u rms L/v — > oo) and if in that limit the 
velocity field u v — > u, for a rough (non-smooth, singu- 
lar) u. In that case, as k, v — > 0, 

G2«{x.f,t f \xo,to) -> G u (x f ,tf\xo,to) (13) 

for a nontrivial probability density G u . The Lagrangian 
trajectories can remain random as K,v — > 0! This phe- 
nomenon has been called "spontaneous stochasticity" 7 , 
because of the analogy with spontaneous symmetry- 
breaking in condensed matter physics and quantum field- 
theory where, for example, a ferromagnet may retain 
a non-vanishing magnetization even in the limit of van- 
ishing external magnetic field. It is important to ap- 
preciate, however, that "spontaneous stochasticity" is a 
very different type of randomness than is usual in turbu- 
lence theory, associated to a random ensemble of velocity 



fields. Instead, the randomness in ( 13 ) is for a fixed (non- 



random) ensemble member u. The limiting distribution 
consists of time-histories that are all solutions of the same 
deterministic initial-value problem 



x = u(x,t), x(i )=x . 



(14) 



As is clear from the limiting probability measure is 
in a certain sense the uniform or equal-weight distribu- 
tion over all such solutions. We note in passing that 
Kneser's theorem in the mathematical theory of ODE's 
implies that, whenever there is more than one such solu- 
tion, then there is in fact a continuous infinity of solutions 
(see HartmarP^, section II. 4). 



The above results have been rigorously established^^ 
for some model problems of turbulent advection, primar- 
ily the Kraichnan model of advection by a Gaussian ran- 
dom velocity field with zero mean and covariance 

«(x,t)^(x',0) = [DoSij - S£(x - x')]<5(t - 0- 

The velocity fields are temporal white-noise and spatially 
rough for l v < r < L u , with a Holder exponent < h < 
1, but smooth for r < l v . We discuss the model only for 
incompressible flow, in which case 



WAT) 



Di[(l + h)Sij - hfifj] 
D 1 l 2 u h - 2 [26 ij - hfj}r 2 



.21 1 



4 « r < L u , 
r < t v . 



(15) 

The velocity realizations u u (x,t) are divergence- free and 
Holder continuous with exponent h for v = D\l 2h — > 0. 
The kinetic energy spectra are power-laws k~ n for the 
"inertial range" 1/L U <C k <C 1/4, with n = 1 + 2h 
and thus 1 < n < 3. The key feature which makes this 
model analytically tractable is the Markovian property 
in time, which leads to the exact validity of Richardson's 
2-particle diffusion equation in the form 



d t P(r,t) 



d_ 



dP 

S^v) — (r,i) ) -! 2,v. .,/>(!•./). 



with r = xi — X2 . Note that Richardson's original equa- 
tion is obtained for h = 2/3 rather than the Kolmogorov 
value h — 1/3, a peculiarity of the white-noise approxi- 
mation. The main features of the solutions of this equa- 
tion can be inferred from a study of the dispersion. Mul- 
tiplying ( 16 ) by r^ri and integrating over r leads to 



j t {r k (t)r e (t)) 



2(S^(v(t)))+iKS u 



(17) 



An analysis of fll6| and (17) leads to the following key 



results for the ensemble of particles which are all started 
at the same point (ro = 0) . 

For Pr = - < 1, the dispersion is that of two indepen- 
dent Brownian motions in three space dimensions 



l/h 



(r 2 (t)) ~ 12nt (18) 

for short times t <C t K — ~ with £ 2 — (jj^j 
£ 2 /(Pr) l ^ h , but of Richardson-type 

(r 2 {t)) ~ g^t) 1 '^ (19) 

for longer times t K <C t <C th = L 2 ^ h ^ /Di, with <?/,. a 
constant independent of k. For Pr > 1, the behavior is 
a bit more complex. The result (JTsJ) still holds at very 
short times t <C t v = l^/v, but at large Prandtl numbers 
there is an intermediate range of exponential growth 



(r 2 (t)) (X nt v e 



2A„t 



(20) 



for t v <C £ <C ln(Pr)i„. Here A„ cx t v 1 is the leading Lya- 
punov exponent of the smooth advecting velocity field. 
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The dispersion at longer times again follows the Richard- 
son law (19). In either case, the Richardson law is valid 



once (r 2 (Fj}> max{^,^} and thus holds for arbitrarily 
small times t > in the limit as v 1 k — > 0. 

A non-vanishing dispersion implies that the La- 
grangian trajectories must stay random in the limit. Al- 
though the diffusion equation ( 16 ) has been averaged 
over velocity realizations u, it must be the case that 
P u (r,i|0, 0) ^ S 3 (r) for t > and for a set of u with 
nonzero probability, or otherwise the average over u 
would also be a delta-function! The physical mechanism 
of spontaneous stochasticity is clearly the "forgetting" of 
the length-scales £ K , i v by Richardson diffusion for suffi- 
ciently long times, t ^> t K ,t v , with those times also van- 
ishing in the limit v, k — > 0. In the case of incompressible 
flow, the limiting distribution is completely independent 
of how the limit is taken. We note in passing that this 
is not true in general for compressible flows and that the 
possibility for Lagrangian particles to stick as well as to 
stochastically split allows there to be different limits, de- 
pending upon the Prandtl number Pr, in the limit as 
v, k —> C) 6 I 10 I 12 I. However, for incompressible flow the lim- 
iting distribution is very universal and robust. 

The same limit is obtained for incompressible flow even 
with ft — or Pr = oo, if the randomness is intro- 
duced through the initial conditions rather than stochas- 
tic noise. Consider the solution of the initial-value prob- 
lem for the Kraichnan ensemble of velocities 

j^E(t) = u"(x(i),t), x(t )=x + ep 

where J) is a zero-mean, unit-variance random vector 
with probability density Q(p). One may interpret e as 
the size of error in measuring the initial position of the 
particle. This corresponds to solving the Richardson dif- 
fusion equation ( |16[ ) with k — and initial condition 
P e (r,t ) = e- 3 (Q * Q)(r/e) so that F £ (r,i ) -> <5 3 (r) as 
e — > 0. However, if the limits are taken v — > first and 
e — > subsequently, then the solution of the Richard- 
son equation does not degenerate to a delta-function for 
t > 0. This may be seen by solving for the dispersion 
from eq. Jl7} with k = and (r 2 (0)) = 2e 2 . When e < l v 



(At)) 



2e 2 e 2X„t 



(21) 



for times t -C t v \n(£ v /e) but for longer times follows the 
Richardson law (19). When l v < e instead, then the 



any time t > if v, e — > with e vanishing slower than 

^ e -o(iv^< 1 - h >)_ 

As this last discussion should make clear, the phe- 
nomenon of "spontaneous stochasticity" is not especially 
connected with the random perturbation of the motion 
equations in (8j). Instead it is the advection of the parti- 
cle by a rough velocity field and the "forgetting" of the 
initial separations which makes the Lagrangian particle 
motions intrinsically stochastic. Spontaneous stochastic- 
ity should not be confused with "chaos" , as that term is 
used in dynamical systems theory^l. For chaotic dynam- 
ical systems with a smooth velocity field, one sees only 
exponential growth of deviations as in eq.|2T|). Because 
this result is proportional to e 2 , the initial separation is 
never "forgotten" and, for all times, (r 2 (t)} — > as e — » 0. 
For chaotic dynamics, any imprecision in the initial data 
is exponentially magnified, leading to loss of predictabil- 
ity at long enough times. Spontaneous stochasticity cor- 
responds instead to \ v = +oo. The solution is unpre- 
dictable for all future times, even with infinitely precise 
knowledge of the initial conditions! 

We remark finally that spontaneous stochasticity has 
very important implications for turbulent advection of a 
passive scalar^. Note that 



d 
dt 



J d 3 x 9 2 (x,t) = -2k J d 3 x |V6>(x,i)| 2 , 



so that, naively, the integral is conserved for k = 0. 
However, in the infinite Reynolds-number, fixed Prandtl- 
number limit {v, k — ¥ with Pr constant) 

0(x,*)= / d 3 x 0(x ,t )G u (x ,t \x,t), 



using the results ( 11 1 and ( 13 ). The scalar at the present 



time is a nontrivial average along stochastic Lagrangian 
trajectories of its values at an earlier time, with molecular 
mixing replaced by turbulent mixing. In particular, the 
scalar intensity — for every velocity realization, without 
ensemble averaging — is still dissipated 



d 6 x9 z {yL,t)< / d A x 6H(x,t ), t > t 



even as v 1 k — > 0! This is the scalar analogue of the 
dissipative anomaly of Onsager for fluid turbulence^HSH. 
The Lagrangian mechanism is spontaneous stochasticity. 



short-time behavior is diffusive 



(r 2 (t)) - 2e 2 + (const.) D^t 



(22) 



C. Experimental and Numerical Results 



for times t 



Dispersion for such a "ran- 



dom clo ud" o f initial positions was first considered by 
BatcheloJ 30 * 31 !, who obtained instead ballistic growth 
cx t 2 for hydrodynamic turbulence. The diffusive re- 
sult above is an artefact of the white-in-ti me v elocity. 



At times t £ ^ - the Richardson law (|19|) holds. 



As in the previous cases, the Richardson law holds for 



Our quantitative discussion above has been based upon 
the original Richardson theory and, in particular, his 
diffusion equation (16). This equation is exact for the 



Kraichnan whitc-in-time velocity ensemble but is only an 
approximation for hydrodynamic turbulence, where sev- 
eral of its quantitative predictions are known to be incor- 
rect. We have already mentioned the diffusive short-time 
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growth in dispersion for a particle cloud, eq.(22|, which 
is ballistic for fluid turbulence. (Note that the ballis- 
tic regime is correctly predicted by Richardson's diffu- 
sion equatio n for a suitably time-dependent ed dy dif- 
fusion tensoJ^^.) The diffusion equation ( 16 ) holds 



in the Kraichnan model also for backward-in-time dis- 
persion. However, actual dispersion rates are different 
forward and backward in time because of the negative 
skewness of turbulent velocity increment s] 35 * 36 !. There is 
presently no quantitative theory of turbulent dispersion 
which successfully accounts for all aspects of the phe- 
nomenon, ft is necessary to stress that the prediction of 
"spontaneous stochasticity" has more general grounds in 
the mathematical theory of ODE's and is not dependent 
upon the diffusion approximation (16). Nevertheless, in 



the absence of any fully successful, quantitative theory, 
it is important to develop understanding from numerical 
simulations and laboratory experiments. We shall here 
briefly review the empirical studies of turbulent disper- 
sion and the status of Richardson's theory. In particular, 
we shall present some new numerical results of our own 
on stochastic particle advection according to eq.® for a 
turbulent velocity field. 

We confine our discussion to just some of the latest 
studies by experiments^® and simulations^^ at the 
highest Reynolds numbers. For more complete surveys 
of the litera ture, see those papers and also recent review 
articled 23 - 12 ^. Athough the t 3 -law (|4) and the stretched- 
exponential PDF (3]) are probably the most famous pre- 
dictions of Richardson's theory, even more important for 
our discussion is the "forgetting" of initial separations. 
If ?*o is the initial particle separation distance and e is 
energy dissipation per mass, then for times much greater 
than t = (rl/e) 1 / 3 both (r 2 (t)) and P(r,t) should be- 
come independent of r$. As we have seen, this is the cru- 
cial physical mechanism underlying spontaneous stochas- 
ticity. In general, it has proved rather difficult to observe 
in a completely consistent and convincing way all of these 
predictions of Richardson's theory. 

Experiments of Ott and Mann at maximum Taylor- 
scale Reynolds number Rer = 107 observed both a i 3 law 
and the Richardson PDF, but varied ro only by a factor 
of 1.5 around the value ro = 10£„ (for l v — (^ 3 /e) 1 / 4 
the Kolmogorov dissipation length). Thus, they provide 
no information on collapse independent of ro. A series 
of experiments by Bodenschatz and collaborators^"^! at 
substantially higher Reynolds numbers up to Rex = 815 
fail to see a t 3 law and instead produce results consis- 
tent with Batchelor's ballistic t 2 range. However, their 
smallest achievable value of the initial separation was 
only about ro = 3Q£ U , so that it is arguable that they 
need longer times and still higher Reynolds numbers. At 
their smallest value of r they did observe Richardson's 
PDF ^ and for r = 20-150£„ they see results roughly 
consistent with Richardson's predictions for the quantity 



,2/3 



(*)) 



2/3 



and also some tendency to collapse inde- 



pendent of ro (see Xu et alP^, Fig. 1). 

Numerical simulations of Ishihara and Kaneda^ at 



R&T = 283 claimed to observe an inertial-range t 3 law 
for ro in a range between 5-45£„, but no tendency what- 
soever for collapse at long times independent of r . How- 
ever, Biferale et alP 2 "! in a simulation at nearly identical 
Rer — 284 and with ro in a range of 1.2-19.6^„ see ex- 
actly the opposite: only a slight indication of a t 3 -law 
for their smallest separation but a strong tendency to 
collapse at long times. They also observe Richardson's 
PDF (|) for r = 1.2£ v . More recently, Sawford et alP 
have performed simulations with maximum Rex = 650. 
They find the best evidence yet of Richardson's predic- 
tions for the dispersion (see their Fig. 4), with a reason- 
able t 3 range for ro = 4J„. Other values of ro in the 
range of 0.25-256^„ do not give a convincing i 3 -law but 
do verify a tendency for collapse at long times. 

In view of the incomplete verification of Richardson's 
predictions, we decided to undertake our own numeri- 
cal investigation. Unlike previous works, however, which 
have nearly all studied deterministic fluid particles with 
variable initial separations ro, we have instead studied 
the problem of stochastic particle advection according to 
the eq.Q. The velocity field is obtained from a 1024 3 
pseudospectral numerical simulation of forced, statisti- 
cally stationary turbulence at Rer — 433. The flow data 
is available online at http://turbulence.pha.jhu.edu 
and fully documented there and in papers^^. The 
entire flow history for about one large-eddy turnover 
time L u /u' is archived at a time resolution suitable for 
particle-tracking experiments, with spatial and tempo- 
ral interpolation implemented within the database. This 
is very convenient for our purposes, since it permits us 
to study particle dispersion backward in time as well as 
forward. As we have discussed above, it is backward dis- 
persion that is most relevant for turbulent mixing^!. 

We have studied two values of the Prandtl number 
Pr = v/k = 1 and 0.1. We solved ^ using the sim- 
plest Euler-Maruyama scheme and also, for convergence 
analysis, an explicit, 1.5th-order strong scheme (Kloeden 
and PlaterfHl Section 11.2). We took time discretiza- 
tion dt — 10~ 3 , which guaranteed that particles moved a 
fraction of l v under both turbulent advection and Brow- 
nian diffusion at each time-step. The velocity field be- 
tween stored data-points was interpolated by 6th-order 
Lagrange polynomials in space and piecewise-cubic Her- 
mite polynomials in time. The results were verified to 
be converged in dt both by comparison with the higher- 
order method and with the Euler scheme at a halved step 
size. For each initial location Xo we evolved N = 1024 
independent particle realizations, giving 523,776 particle 
pairs, over the whole time range of the database. For 
forward tracking we started particles at to = and for 
backward tracking at tf = 2.048 (the final time in the 
database). We then averaged all results over 256 initial 
locations Xo obtained by choosing 4 independent, uni- 
formly distributed points from each of 4 3 = 64 subcubes 
of the whole flow domain. 

Our results for particle dispersion are given in Fig.l. 
We present there a log-log plot of (r 2 (t)) (normalized by 
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(a) Backward Dispersion 



- Pr = 1 
-12* 
■1.35i 3 
■Pr = 0.1- 
■120t h 




time t viscous units 



(b) Forward Dispersion 




time t [viscous units] 



FIG. 1. Mean dispersion of particle pairs: (a) backward dis- 
persion, (b) forward dispersion. The backward in time results 
are plotted against t' = tf — t and all quantities are non- 
dimensionalized with viscous units (see text). The Pr = 1 

results are plotted with solid lines ( ), Pr = 0.1 results 

with dashed lines (— — — ). We color code the lines with 
blue for raw data, green for short-time molecular diffusion, 
and red for long-time Richardson diffusion. 



£ 2 ) versus time t (normalized by t v = £l/v). For back- 
ward dispersion we take t — > t' = tf — t, to facilitate 
comparison with the forward-in-time results. In the vis- 
cous uni ts t hat we employ, the early-time diffusive sep- 
aration (181 becomes (r 2 (t)) ~ YitjPr. This regime is 



clearly seen for both backward and forward cases and 
for both Pr = 1 and 0.1. Furthermore, for the unit 
Prandtl number cases, we see a convincing transition to 



a i 3 -law for t > 1. This occurs slightly earlier for back- 
ward dispersion than for forward. Also, we find that 
the asymptotic Richardson-Obukhov constant is greater 
for backward dispersion than for forward, in agreement 
with earlier results^. An average of the local constants 
g(t) = d(r 2 (t)) / d(t 3 ) in the i 3 -scaling range gives val- 
ues of 50 = 1-35 for backward dispersion and go = 0.64 
for forward dispersion. The latter agrees perfectly with 
a recent theoretical prediction^ anc j both are generally 
consistent with previous values^. Pure cube scaling laws 
with these coefficients are plotted in Fig.l for compari- 
son with the numerically obtained mean dispersions. The 
agreement is obviously quite good at long times, espe- 
cially for the backward case. Even more importantly, 
the Pr = 0.1 dispersions show a very clear trend to ap- 
proach the same cubic laws at sufficiently large times 
t 3> t K — (k 3 /e) 1 / 4 or, in viscous units, t ^> (Pr) -3 / 4 . 
Our Fig.l thus provides strong evidence of the "forget- 
ting" of the molecular diffusion time-scale t K by turbulent 
Richardson diffusion at long times, which is the essential 
ingredient of spontaneous stochasticity. 

To be completely conclusive, we would need to see col- 
lapse of the dispersion curves for different Pr in a range 
where both show t 3 -scaling. We see no clear Richard- 
son t 3 for the Pr = 0.1 cases in the time ranges plot- 
ted. We cannot continue the time-integration further for 
two reasons. First, the velocity field from the turbu- 
lence archive contains no data for longer times. Second, 
the rms dispersion distance L(t) — J (r 2 (t)) at the fi- 
nal time has reached a value L(tf) = 1, just slightly 
smaller than the velocity integral length-scale L u = 1.376 
for the flow. (This is expected, since tf is about one 
large-eddy turnover time.) To integrate further to see a 
conclusive collapse, we would need a numerical simula- 
tion at higher Reynolds numbers and integrated to longer 
times. However, our confirmation of Richardson diffusion 
at Rer = 433 using stochastic Lagrangian trajectories is 
comparable to, or even better than, the results of Saw- 
ford et alPat i?e T = 650 using deterministic Lagrangian 
trajectories. (See Fig. 4 in that paper). There are two 
plausible arguments why this should be so. In the first 
place, using stochastic trajectories, all the particles start 
at the same point. At the time t « t u when (r 2 (f) « £ 2 
(for Pr = 1), the particles are not randomly placed in the 
flow, but have already been experiencing relative advec- 
tion by different-sized eddies at the onset of the inertial 
range. Thus, they begin to experience Richardson dif- 
fusion at that time. However, using the usual technique 
of seeding the flow with particles at initial separations 
r o = C one would still need to wait some additional 
time for the initial configuration to be "forgotten". A 
second reason is that backward dispersion is faster than 
forward dispersion, so that the range of t 3 scaling occurs 
even earlier in that case. The technique of stochastic 
Lagrangian trajectories appears to be promising in the 
numerical study of Richardson diffusion. 

In order to make a completely convincing case that 
we are observing Richardson diffusion, we have also nu- 
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merically calculated the PDF P(r, t) of the particle- 
separations. Our results for Pr = 1 are presented in 
Fig. 2, with the normalizatio n dr r 2 P{r,t) — 1. As 
has been previously observe d 23 * 37 * 42 ! , Richardson's ana- 
lytical formula for the long-time PDF of separation dis- 
tances implies that all the PDF's at differ ent times will 
collapse when scaled with L(t) = J (r 2 (t)). In fact, equa- 
tion ([3| is equivalent to 

2/3 



Backward Dispersion 



L 3 (t)P(r,t) = exp 



L(t) 



with numerical values 

a = (1287/8) 1/3 = 5.4387 

and 



4.7617. 



Thus, Richardson's theory makes a parameter-free pre- 
diction that a log-linear plot of L 3 (t)P(r,t) versus 
(r/L(t)) 2 / 3 should give a straight line with slope —a and 
j/-intercept /3. In Fig. 2, therefore, we have plotted our 
PDF's in this way, at three times t = 22.37, 33.57, 44.79 
all lying in the range of t 3 scaling. We have also plot- 
ted the straight line predicted by Richardson's theory. 
We see that the PDF's scaled in this way collapse very 
nicely. Furthermore, except for some deviation at small 
r in the backward dispersion case, they very closely agree 
with the predictions of Richardson's theory. 



III. STOCHASTIC FLUX-FREEZING 

The standard views on flux-freezing in high- 
conductivity plasmas are inconsistent with the phe- 
nomenon of spontaneous stochasticity. It is nearly ubiq- 
uitously argued that flux-freezing should hold better as 
magnetic diffusivity A — > 0. However, high magnetic 
Reynolds numbers are usually associated also with high 
kinetic Reynolds numbers. If kinematic viscosity v — > 
simultaneously with the resistivity and the plasma be- 
comes turbulent, then Lagrangian trajectories will no 
longer be unique. Which fluid trajectory shall a mag- 
netic field line follow if there are infinitely many such 
trajectories? This is the paradox of flux-freezing. 

As we shall argue below, a form of flux-freezing does 
survive at small resistivities and viscosities, but in a novel 
stochastic sense. Before we make this argument, however, 
we shall first discuss the related subject of flux-freezing 
properties of resistive hydromagnetics. 



A. Resistive Hydromagnetics 

In this subsection we shall discuss magnetic fields that 
satisfy the resistive induction equation 




(r/L{t)fl 3 



Forward Dispersion 




d t B = Vx(uxB - AVxB), 



(r/L(t)) 2 / 3 



FIG. 2. Probability densities of pair-separation distances: 
(a) backward dispersion, (b) forward dispersion. The quanti- 
ties at different times are normalized by Lit) = \J (r 2 (t)) as 
shown. We plot densities for three different times in the t 3 
scaling range, with blue for t = 22.37, green for t = 33.57, 
and red for t = 44.79. In the backward case these times cor- 
respond to t' =%] — %. The solid black line gives Richardson's 
analytical prediction for the density (see text). 



with A = rjc/^TT the magnetic diffusivity. It is impor- 
tant to stress that our analysis applies here applies to a 
very general velocity field u. It may be incompressible or 
compressible. It may be externally prescribed or it may 
satisfy a dynamical equation that contains B itself. For 
example, u may be the plasma velocity that obeys the 
standard magnetohydrodynamic momentum equation or 

— VxB, the electron 



(23) it may be taken to be u e = u 
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fluid velocity in Hall magnetohydrodynamiciP^. Our only 
assumption in this section shall be that u is a smooth 
vector field. 

A priori there is no obvious way how to describe 
magnetic-line motion for non-ideal plasmas. One ap- 
proach that has been widely employed in discussions of 
magnetic reconnectiorP^ (or Kuslrud^, Section 3.4) is 
to introduce a "slip velocity" Au = Bx(AVxB)/5 2 . 
In that case, one may attempt to introduce an "effec- 
tive velocity" u* = u + Au of the field lines. Unfortu- 
nately, this approach is not generally successful because 
(Au)xB = -AVxB if and only if B-E = (or B J = 
for Ohmic non-ideality). As has been emphasize d 48 * 49 *, 
no effective-velocity approach is satisfactory for discus- 
sions of three-dimensional magnetic reconnection. In 
fact, those authors show that, even if the non-ideality is 
spatially localized, there generally exists no smooth ve- 
locity field whatsoever such that <9 t B = Vx(u t xB) 
for a non-ideal plasma. 



For magnetic fields that obey (23), however, there is 



a natural and consistent way to describe line-motion as 
a process of stochastic advection. Such approaches have 
already been employed for some time in discussion of 
kinematic ma gneti c dynamos, at least for incompressible 
velocity field d 50 * 5 M Recently, we gave a rigorous proof 
of stochastic flux-conservation properties for nonlinear 
hydromagnetic models using mathematical methods of 
stochastic analysis^!. We shall present here a more phys- 
ical demonstration of these results using path-integral 
methods which, also, extends their validity to compress- 
ible fluid models. 



To begin, we note that the induction equation ( 23 1 may 
be rewritten as 

diB + (u-V)B = (B-V)u - B(V-u) + AAB. (24) 

In this form it is the same as the scalar advection equa- 
tion (10), except for the additional two terms on the 
righthand side. The path- integral formula (11) for the 



scalar solution may thus be easily adapted to this situa- 
tion. The solution of (24) with initial condition B(to) = 
B is given by the "sum-over-histories" formula 



B(x,t) = f 

J a 



Vsl B (a(t ))-J"(a,t) 



eXp ("4l/ ^|a(r)-u>(T),T)| 2 ) 



(25) 



where jT(a, r) where is a 3 x 3 matrix and B is interpreted 
as a 3-dimensional row vector. J~ satisfies the following 
ODE along the trajectory a(r): 



d_ 



J(n, t) = J(sl, r)V B u(a(r),T) 



^(a,T)(V x .u)(a(r),r), (26) 



with initial condition J7"(a, to) = I. It is easy to check 
by taking the time-derivative of (25) and using (26) to 



show that the induction eq.(24) is satisfied. Just as for 



the scalar problem, the condition a(t) = x on the path- 
integral trajectories implies that they correspond to so- 
lutions of the stochastic equation 



(It 



a(r) = u(a,r) + V2At7(t), a(i) = x 



(27) 



integrated backward in time from t = t to r = to. 

However, the stochastic equation (27) may also be in- 
tegrated forward in time from r = to to r = t. In that 
case, the same ensemble of trajectories may be obtained 
by considering only those particles with initial locations 
carefully selected to arrive at x at time t, for a given 
realization of the white-noise rj(t). With a slight change 
of notation, we may characterize this ensemble of time- 
histories as those x(t) which solve 



|:x(a,T) = u(x(a,r),r) + V2At7(t), t > t ^ 
x(a,t ) = a 



such that the inverse map a(x, t) to x(a, r) specifies 
the starting point by a = a(x,t). Notice that (28) is a 



stochastic generalization of the usual equation for a La- 
grangian flow map x(a, t) of a particle with initial "label" 
a and that a(x, t) is the "back-to-labels" map. It is easy 
to show furthermore by applying V Q to (28) that 



J(a,t) = \ V a x(a,t). 

det (V a x(a,t)) 



(29) 



solves equation (26) with initial condition J7"(a, t ) = I. 



It is therefore possible to re-express the path-integral for- 
mula (25) as 



B(x,t) = 



1 



det(V a x(a,t)) 



B (a)-V a x(a,t)| 5(X]t) 



_ ( 30) 
The overbar • represents the average over realizations 

of the random white noise process rj(t) in (28). 



We shall call the above result the stochastic Lundquist 
formula, since it is the stochastic generalization of the 
standard Lundquist formula^ (or KulsrucP, section 4.8). 
It may be cast into a more familiar form by noting that 
the determinant that appears there can be interpreted as 
the ratio of initial and final mass densitieJ^ 



det(V„x(a,t)) 



Po(a) 



p(x(a, t),t) 



It follows that the vector field B/p is stochastically 
"frozen-in" and advected along stochastic Lagrangian 
trajectories, where B is defined to be the quantity under 
the overbar in (30). Notice, therefore, that the average 



in (30) is not over the frozen-in field B/p, but rather over 



the magnetic field B itself. This is necessary in order to 
reproduce the Laplacian term in ( 24 ) , which has the form 
AAB and not A A (B/p). 

The use of the stochastic Lunquist formula is illus- 
trated in Fig. 3. The aim is to calculate the magnetic field 
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average 




FIG. 3. Illustration of the stochastic Lundquist formula. 
Three stochastic Lagrangian trajectories running backward 
in time from a common point are shown in red, green, blue. 
Starting field vectors, represented by correspondingly colored 
arrows, are transported along the trajectories, stretched and 
rotated, to the common final point. These are then averaged 
to give the resultant magnetic field at that point, indicated 
by the black arrow. 



B at spacetime point (x, t) . The first step is to generate 
an ensemble of stochastic Lagrangian trajectories solving 
(27) backward in time from x at time t to random loca- 



tions a(to ) at the initial time to- The path-integral for- 
mula (25) sums over all such random time-histories. We 



show in Fig. 3 (top) three stochastic trajectories gener- 
ated numerically from the turbulence database together 
with the starting magnetic field vectors B , indicated by 
arrows at the starting locations a(io). The next step is to 
transport each of the field vectors in the usual "frozen- 
in" fashion along the stochastic Lagrangian trajectories 
to the final spacetime point (x, t) . The result is an ensem- 
ble of field vectors B at that point, stretched and rotated 
by the flow. These are illustrated in Fig. 3 (middle) by 
the collection of three arrows at (x, t), obtained by trans- 



porting the three initial vectors. In the usual determinis- 
tic Lundquist formula there would be just one trajectory 
and one vector B at the final point, which would give the 
desired magnetic field. Now however as the final step one 
must average over the ensemble of random vectors B in 
order to obtain the resultant magnetic field B(x, t). This 
is illustrated by the black arrow in Fig. 3 (bottom). In 
contrast to the previous transport step which preserved 
line-topology (in each individual realization), the final 
averaging step resistively "glues" the transported lines 
together and changes the magnetic field-line topology. 

There is an elegant reformulation of the stochastic 
Lundquist formula which must be mentioned here, both 
because of its conceptual simplicity and also because 
of its potentially greater generality (see next subsec- 
tion). Consider any smooth, oriented surface S at fi- 
nal time t. Then the formula (30) may integrated in x 



over the surface S, with respect to the vector area el- 
ement dA(x) = dxxdx, and the ensemble-average and 
surface-integration interchanged on the righthand side. 
Because the expression under the overbar is the one that 
appears in the usual Lundqust formula, the standard 
multi-variable calculus manipulations convert this into 
a surface integral over a.(S,t), the surface S randomly 
advected backward in time to the initial time to. As 
before a(-,t) = x _1 (-,t) is the "back-to-label map" for 
the stochastic forward flow. The result is the following 
stochastic Alfven theorem 



/ B(x,t)-dA(x) = / B (a)-dA(a) , t > t . 

r-i (31) 

This result generalizes a previous theoreirPU to compress- 
ible plasmas. Eq.(31 1 expresses the conservation of mag- 
netic flux on average, as illustrated in Fig. 4. An ini- 
tial loop C, boundary of the surface S, is shown there in 
black. This is stochastically advected backward in time 
to give an infinite ensemble of loops at the initial time 
t - These are represented by the three colored loops. The 
ensemble-average of the magnetic flux through the collec- 
tion of loops at the initial time to is equal to the magnetic 
flux through the loop C at the final time t. 

The stochastic Alfven theorem is an example of what is 
called a "martingale property" in probability theory. The 
magnetic flux through each advected loop at the earlier 
time to is unequal to the magnetic flux through C at time 
t. Nevertheless, the mean flux remains the same. Note 
that this result implies an irreversibility or an "arrow of 
time" since it only holds for backward stochastic advec- 
tion of loops. Backward-in-time is the causal direction, 
since the magnetic flux at the present must be obtained 
as an average of past values and not of future values. 
If we assumed a "forward martingale" property then we 
would obtain instead the magnetic induction equation 
(24) with a negative resistivity term — AAB. Note, in 



fact, that the stochastic Alfven theorem (backward in 
time) is mathemati cally eq uiva lent to the usual resistive 
induction equation (23) or (24) 52 . 
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FIG. 4. Illustration of the stochastic Alfven theorem. Shown 
are three members (red, green, blue) of the infinite ensemble 
of loops obtained by stochastic advection of a loop C (black) 
at time t backward in time to to ■ The average of the magnetic 
flux through the ensemble of loops is equal to the magnetic 
flux through C. 



B. High-Reynolds-Number Limit 

We now consider the limit of large kinematic and mag- 
netic Reynolds numbers. For simplicity we shall assume 
that Pr m = v j\ remains fixed as v, A — > 0. 



Consider the Feynman-Kac formula (25 1. By a naive 
application of the Laplace method, one would assume 
that the path-integral collapses to a single determinis- 
tic trajectory as A — > 0, with rms fluctuations of or- 
der (At) 1 / 2 for small but nonzero A. This is precisely 
the heuristic estimate of line-slippage made by Kulsrud 2 
which was quoted in the Introduction. This estimate is 
rigorously correct if the velocity and magnetic fields are 
assumed to remain smooth in the limit v, A —¥ 0. Thus, 
the heuristic estimate is correct if the plasma flow re- 
mains laminar, but this will be the exception rather than 
the rule at high Reynolds numbers. In a turbulent flow, 
the behavior will be quite different. As we can see from 
our Fig. 1 for incompressible hydrodynamic turbulence, 
the heuristic estimate is only valid for very short times 
smaller than the resistive time t\ = (A 3 /e) 1 / 4 . At longer 
times, the rms slip distance of the field lines instead fol- 
lows the Richardson law ~ (et 3 ) 1 / 2 , independent of v 
and A. The quantitative behavior will be different in plas- 
mas with strong magnetic fields, due to the effects of the 
Lorentz force, as discussed more in section V. However, 
the qualitative behavior must be the same whenever the 
advecting velocity is turbulent and spatially rough. 



The Feynman-Kac formula (25) is not very well-suited 
to analyzing the limit of high Reynolds numbers, how- 
ever, because the velocity-gradients that appear in the 
definition of the matrix diverge in that limit. Like- 
wise, the gradients of the Lagrangian flow map that ap- 
pear in the definition (29) of J" are expected to diverge. 



Alfven theorem (31 ), is more likely to remain meaningful 
in the limit of infinite Reynolds number. The backward- 
advected loops a A (C, t) at finite values of A, v are ex- 
pected to approach well-defined curves a(C, t) as A — > 0, 
which, however, are not rectifiable but fractaPSHSfll To 
make mathematical sense of magnetic flux through such 
fractal loops, we may introduce the vector potential 
Ao = (curl) _1 Bo and rewrite the flux through surface 
a(S, t) as a line-integral around its perimeter a(C, t) . We 
may then further transform by change of variables to a 
line-integral around the original loop C, as: 

d> A (a)-e?a = d> Ao(a(x, i))-<fa(x, t). 

Ja{C,t) JC 

The integral on the right may be interpreted as a gener- 
alized Stieltjes integral, which is well-define d as l ong as 
the map a(x, t) is suitably Holder continuou d 61162 '. 

It may seem from our arguments to this point that the 
validity at very high Reynolds numbers of the stochastic 
flux-freezing result (31 ) is dependent upon the particular 



stochastic representation of resistive effects employed in 
eqs.(28) and (30). However, the same "martingale prop- 
erty" can be obtained in the limit v, A — > by a differ- 
ent arg ument that employs only the standard Lagrangian 
flow 63 -'. Define the deterministic flow, as usual, by 



d x^a,T)-u-(x"(a,T),r) 



dr 



T>t n 



x"(a,i ) = a 



(32) 



where the superscript v is a reminder that the dynam- 
ical equation for the advecting velocity field contains a 
certain viscosity v — (Pr)\. Correspondingly one defines 
the inverse map a" (;t) = (x")" 1 ^,*). Stochasticity can 
be introduced by assuming small random perturbations 
of the loop, taking C — > C + eC, where C is a ran- 
dom loop from a well-behaved ensemble^. Thus, e can 
be regarded as the spatial-resolution in determining the 
precise form of the loop C. We may then argue that, 
at least for incompressible flow, the ensemble of loops 
a A (C, t) obtained from stochastic advection in the limit 
A — > coincides with the ensemble of loops obtained 
from deterministic advection a v (C + eC,t) taking the 
limits first v 1 A — > and then e — > 0. As discussed in 
the previous section, this is rigorously known to be true 
for point-particles advected by veloc ities selected from 
the Kraichnan white-in-time ensembl e! 9 ' 12 '. The physical 
mechanism is just turbulent Richardson diffusion. We 
therefore conjecture that the same result holds for loops. 
If this is so, then the double limit 



lim lim 

e->0 jaA->0 



'(C+eC,t) 



Ao(a)-<ia 



gives precisely the same ensemble of fluxes with the same 
distribution as in the previous approach. In that case, 
(31 1 must again hold in the limit, or, more precisely, 



The integrated form of flux-conservation, the stochastic 



f A(x, £)-dx = lim lim (p A (a)-da, 



(33) 
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where A = lim^^-^o A A and the overbar now indicates 
average over the ensemble of loop perturbations C . This 
is a nontrivial result because for an individual loop C 



A. Lagrangian Description of Dynamo 



d 
df 



f A A (x,t)-dx = -A f (VxB A )(x,t)-dx 



Since VxB A diverges in the limit i/,X — > 0, there is no 
reason to expect that the righthand size vanishes in that 
limit. This is the standard argument how flux- freezing 
can be violated in thin current-sheets. It stands to rea- 
son that magnetic flux through an individual Lagrangian 
loop will fluctuate in time and not be conserved. Nev- 
ertheless, our arguments lead us to conclude that flux- 
freezing in turbulent flow is still preserved in the mean 
sense (33) at infinite Reynolds number. 



A scalar resistivity of the form in (23) or (24), in fact, 



plays no essential role in our arguments. Any microscopic 
plasma mechanism of "line slippage" will be accelerated 
by turbulent advection as soon as the lines have sepa- 
rated by a distance of order i v , the viscous length. More 
realistic mechanisms of line-slippage such as anisotropic 
resistivity in the Braginski equations^ (also KulsrucP, 
Ch.8) or the Hall effect often invoked in theories of fast 
reconnectiorP^ may all serve the same role. After a very 
short time the microscopic plasma mechanism of line- 
slippage, whatever it may be, will be "forgotten" and 
replaced by turbulent Richardson diffusion. 



IV. TURBULENT MAGNETIC DYNAMO 

The stochasticity of flux-freezing plays an essential 
role in the operation of the turbulent magnetic dynamo. 
We have already made a detailed analysis of this in the 
Kazantsev-Kraichnan model of fluctuation kinematic dy- 
namo for a non-helical, inc ompre ssible velocity field, with 
Re m = oo and Pr m = (PGU. Here we shall present 
a more general study. We first discuss how our path- 
integral approach relates to th e standa rd Lagrangian for- 
mulations of magnetic dynamJ 50 * 51 * 67 * ^, in a framework 
that encompasses helical and/or compressible flows and 
turbulent velocity fields with realistic time-correlations. 
We focus on the kinematic dynamo here but much of 
our discussion carries over also to the nonlinear dy- 
namo (which is considered more specifically in section 
V). We then present numerical results for a particular 
case, the kinematic fluctuation dynamo in a non-helical, 
incompressible turbulent velocity field for Pr m = 1, us- 
ing the same hydrodynamic turbulence database that 
was employed in our study of Richardson diffusion (sec- 
tion II. C). We com pare our results with earlier numeri- 
cal studies^ 21 * 72 * 73 ! at lower Reynolds number and, also, 
with our previous analytical results in the Kazantsev- 
Kraichnan model at Pr m — C 15 . As we shall see, sponta- 
neous stochasticity and Richardson diffusion play a very 
similar role in the fluctuation dynamo for both Pr m = 
and Pr m = 1. 



The Feynman-Kac formula ( 25 ) for the magnetic field 
may be rewritten as 

B\K,t) = /' ( i 3 aS fe (a)^(a,to|x,t;u) (34) 

with the definition 



^(Mo|x,t;u)= / Va5 6 {a{t Q )-a)J^cx,t) 
Ja(t)=x 



)r ■ (35) 



This latter quantity is a generalization to a compressible 
flow of the (Eulerian) magnetic Green's function consid- 
ered by Lerche^ and KraichnarP^, expressed as a La- 
grangian path-integral. It completely encodes all the ef- 
fects of the advecting flow. 

A description of the mean-field dynamo is obtained if 
one averages over the ensemble of velocity fields and the 
random initial conditions of the magnetic field. Assuming 
that these are statistically independent (which requires 
that the effects of the Lorenz force be negligible), 



(S*(x,t)> = J d 3 a (B fe (a))^,(a,i |x,0 



(36) 



with Fl(&,to\x,t) = (Fl(a,t \x,t)) the mean magnetic 
Green's function. The same result holds without the 
assumption of kinematic dynamo if the initial magnetic 
field is non-random and the mean Green's function is de- 
fined by a conditional average for fixed BcpM Of course, 
in that case the mean Green's function becomes depen- 
dent upon the magnetic field. The mean Green's function 
involves Taylor 1-particle diffusion, or absolute diffusion 
with respect to a starting point x at time t, with the 
stochastic fluid particle moving backward to time to. In 
the case of homogeneous velocity statistics, 

Fl(p;t ,t) = Fi(p,t \0,t) 

becomes a function of the single variable p = a — x. The 
large- ;0 behavior of the mean magnetic Green's function is 
well-known if the velocity statistics are also isotropic (but 
reflection non-symmetic), in which case the usual a and j3 
effects of mean-field electrodynamics determine the large- 
distance decay (see KraichnarP^, eq.(3.17)). In particu- 
lar, the mean Green's function is non- negligible only for 
p = O ({f3\t - i |) 1/2 ) , with the eddy-diffusivity of the 
mean magnetic field. If the magnetic field statistics are 
also homogeneous, then the mean field (B(x, t)) = (Bo) 
becomes time- independent (no mean-field dynamo), in 
which case ( 36 ) yields the sum-rule 



d 3 P Fl(p;t ,t)=S' k 



(37) 



The exact Lagrangian formulas for the mean-field elec- 
trodynamics coefficients, a,/3, etc., that were derived 
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by Moll'att" 7 and KraichnarPE^ a t infinite conductiv- 
ity, hold also for po sitive re sistivity within the present 
stochastic framework ! 50 ! 51 * 71 ! These formulas involve the 
stochastic displacement field £(x, t) = x — a(x, t) of 1- 
particle turbulent diffusion^]. 

As a side remark, we note that the results on the mean 
Green's function in homogeneous, isotropic turbulence 
which we reviewed above may have limited relevance to 
the description of astrophysical dynamos and laboratory 
dynamo experiments. The separation of scales required 
for the applicability of mean-field electrodynamics often 
does not occur in practice and large-scale magnetic fields 
might not be understood without reference to object- 
specific features, global flow geometry and boundary- 
conditions. In fact, as a general rule, mean-field dynamo 
effect requires not just turbulent diffusion of field-lines 
but also a globally organized motion of magnetic fields. 
To show this, we present an argument based on Faraday's 
law for the mean field 

c> t (B)+V-£ = 0, 

rewritten as a local conservation law for the magnetic 
field vector. Here 



' /!' - </•'/)' + A 



( d{B l ) _ d{Bi_ 
V 9xj 



dxi 



represents the spatial flux of the jth component of mean 
magnetic field in the ith coordinate direction. The space- 
average of (B(x, t)} over a volume V can only change 
in time by a transport of magnetic field lines through 
its surface dV. Furthermore, a simple calculation with 
Ampere's law (J) = ^Vx(B) gives 



d 
dt 



J ±-\(B( x ,t))\ 2 d i x= /'y,,, -r y.'* if'. i'. 



where the integral is over all of space. We see that energy 
in the mean-field (B) grows when the transport of mag- 
netic flux across the closed lines of mean electric current 
reinforces the Amperian fields induced by those currents 
(using the righthand rule). We thus see that mean-field 
dynamo action requires coherent motion of magnetic field 
lines, coordinated over large spatial scales. Note that 
very similar ideas are widely used in condensed matter 
physics to explain, for example, the decay of magnetic 
flux through a superconducting ring by "phase-slippage" 
of quantized field lines^EU. 

Returning to the discussion of small-scale turbulence, 
we note that a formula for the magnetic correlation func- 
tion analogous to ( 36 ) for the mean field can be derived, 



under precisely the same assumptions: 

(5 , (x,t)B 3 (x',t)) =J d3a J rfV (^o(a)^(a')) 

xF«(a,a',i |x,x',t) (38) 

witrP 

F$(a,a',i |x,x',i) = (F*(a, t |x, t)F/(a', t |x, t)). 

(39) 



The behavior of this "two-body" Green's function is 
determined by the properties of turbulent 2-particle 
(Richardson) diffusion effects. E.g. setting x = x' in 
the above formula leads to an expression for mean mag- 
netic energy density (B 2 (x,t)), which is related to pairs 
of stochastic Lagrangian trajectories (with independent 
Brownian motions) that start at points a and a' at time 
to and both end at x at time t. 



Combining ( 36 ) and ( 38 ) gives a formula for the cor- 



relation of the magnetic fluctuations: 
(B\ X ,t)B\ X ',t)) - (^(x,i))(i?V,t)) 

d 3 a J d 3 a' [(^'(a)^(a')> - (B fc (a))(^(a'))] 
xF^(a,a',i |x,x',i) 

+ Jd 3 aj dV(5 fc (a))(^(a')) 

x [F^(a,a',to|x,x / ,t)-F«(a,to|x ) t)F/(a',t |x,t)] .(40) 

The first term on the righthand side represents fluctu- 
ation dynamo due to growth of magnetic fluctuations, 
whereas the second term represents magnetic induction, 
or the generation of magnetic fluctuations from the mean 
field by random advection. Note that for |x — x'| ^> L u , 
the integral correlation length of the velocity field, 

F%, (a, a', t |x, x', t) ~ (a, t |x, t)Ff (a' , t |x, t), 

because the two stochastic particle trajectories become 
statistically independent. As a consequence, the second 



magnetic induction term in ( 40 1 always goes to zero for 
|x — x'| — > oo. The first term will also vanish in that limit 
if <B*(a)i?<(a')> - <B fc (a)>(£<(a')) for |a-a'| oo 
(statistical "clustering" of initial data) . 

The above formulas simplify in the special case of spa- 
tially homogeneous statistics for both the velocity and 
magnetic fields. In particular, (38) becomes 



(B*(r,t)Bi(0,t))= J d 3 p(B h (p)B e Q (0))Fll(p,t \r,t). 

. (41) 

with the homogeneous 2-body mean Green's function 
Fg(p,to\r,t) = J d 3 aF l k j e (a,a + p, t |x, x + r, t). (42) 

For |r| > L u , 

F%(p,to\T,t) ~ J d 3 aFl{a + p;t a ,t)F^n + r;t ,t), 

and the 2-body Green's function is non-negligible only 
for \p-r\=0 ((/3\t - £ |) 1/2 ) • Then implies that 



lim 

|r|— foo 



d 3 pF^(p,t \r,t)=6l5i. 



(43) 



These properties will be used in our discussion of the 
numerical results below. 
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B. Numerical Study of Kinematic Dynamo 

We now present a numerical study of small-scale turbu- 
lent kinematic dynamo at Pr m = 1. We employ the same 
database of non-helical, incompressible fluid turbulence 
that was used in our investigation of Richardson diffusion 
in section II. C. This might appear to be a poor choice at 
first sight, since the conventional view^LLUis ^ e k me _ 
matic fluctuation dynamo at Pr m = 1 is a phenomenon 
of sub- viscous scales. Schekochihin et al. conclude ex- 
plicitly: "This dynamo is driven by the viscous-scale ed- 
dies and whatever the inertial-range velocities might do 
is guaranteed to happen at a slower rate"P For this rea- 
son, the previous numerical studies have taken special 
pains to resolve well the viscous range, at a sacrifice of 
Reynolds number. For example, the highest resolution 
1024 3 simulation of Haugen et alP had a Taylor-scale 
Reynolds number Rer — 230 which is nearly half t hat of 
the database that we employ, for which Rex — 432l 19 | 20 l 
The previous numerical results seemed to verify the idea 
that viscous scales played the dominant role; for example, 
the magnetic energy spectrum in the kinematic regime 
was found to be peaked at wavenumbers a little higher 
than the viscous Kolmogorov wavenumber k v . See Hau- 
gen et alP, Fig. 4 and Schekochihin et alP, Fig. 22(a). 
The viscous range is not so well-resolved in the database 
that we employ, with the grid spacing Ax of the simu- 
lation being slightly greater than 2l v . Nevertheless, our 
study was designed to show the critical role of inertial- 
range advection to the small-scale turbulent dynamo at 
high Reynolds numbers and, thus, was forced to sacrifice 
resolution of the viscous range. Our results will show 
that both ranges play a critical role at Pr m — 1 and, we 
will argue, even for Pr m much larger. 



1. Methods 

Our Lagrangian numerical approach is based upon the 
results in section III. A. We construct an ensemble of 



stochastic particles that solve eq.(27) backward in time 



from common starting point x. The Feynman-Kac for- 
mula (pBl then yields 



where the double-overline indicates an average over two 
ensembles of trajectories a(t), a'(t) with independent re- 
alizations of the Brownian noise. The matrix satisfies 



equation (26) for V^-u = 
d 



dt f 



The exact solution of this equation is an anti-time- 
ordered exponential of the velocity-gradient from t to 
tf and we have here indicated explicitly the dependence 



of J7" upon both times. As a matter of fact, it is numeri- 
cally easier to use the ODE in the initial time to 



dr 



J(a,t f ,T) = -V x u(a(r),T).7(M / ,T) ) J(a,t f ,t f ) = 

which may be solved backward in time from r = tt to 
t = t along with the stochastic equations (27 1. We then 



average over an ensemble of initial conditions Bo (the 
same for each r) to obtain the mean magnetic energy 

(B 2 ( X ,t f )) T = 

— i 

<B*(a(T))Bg(a'(r))) tf, r)J T {a>, t f , r) 

In effect we are solving for the growth of magnetic field 
by moving the time r of the initial conditions backward 
rather than advancing tf forward. Thus, our results be- 
low shall be plotted with respect to the difference variable 
t = tf — t. Assuming ergodicity, an average over space 

(5 2 (i))^i Jd 3 x (B 2 (x,t f )) T 

is equivalent to an average over an ensemble of velocities. 
In a statistical steady state, this average should indeed 
be a function only of the difference variable t = tf — r. 

To further simplify matters, we take as our initial seed 
field for the dynamo a spatially uniform magnetic field Bo 
which is still random, however, and statistically isotropic. 
The covariance choice (BqBq) = |(5 IJ implies a magnetic 
energy initially equal to one. This is not very small, 
but there is no requirement of small field strength in our 
kinematic problem. The formula for the mean magnetic 
energy then factorizes as 



B 2 (x,t f )) T = hr (j( a ,tf,T)J T ( a <,tf,T] 



hr (j{tf,T)J T {t f , 



(45) 



J{a,tf,t ) = t 7(a,t / ,t )V 2; u(a(t / ),i / ), J{a,t ,t ) 



where J(tf,r) = J(a,t f ,T) . Note that F = J J is 
a positive-definite, symmetric matrix which formally re- 
duces in the limit of vanishing noise to the usual (left) 
Cauchy-Green or Finger deformation tensor of continuum 
mechanics. The initially uniform magnetic field does not 
stay uniform but develops small-scale fluctuations by an 
induction effect. There is, in fact, no very precise dis- 
tinction between "magnetic induction" and "fluctuation 
dynamo", as we have discussed elsewhere^, and weak 
uniform seed fields have been used in many previous stud- 
ies of turbulent magnetic dynamo ^ 2 * 73 * 80 !. Thus, at fixed 
r, the magnetic correlation function with this initial seed 
field is dominated at long times t by the leading dynamo 
eigenmode £, 

j (B l {r,t)B j {0,t)) - (const. )e 7t £ y '(r), t ->■ oo 

with 7 the dynamo growth rate. In the opposite limit of 
large distances for fixed t it follows from (41), (43) that 



(B i (r J t)B^(0,i))~^ J ', 



oo 
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for our choice of initial seed field. 

We implemented this scheme numerically by solving 
the SDE p7| backward in time for TV = 1024 samples 
a n (r), n = 1,2, ...N, all started from point x at time 
tf with independent realizations of the noise. We took 
tf = 1.5 and t — 0.5, because the spatially-averaged 
energy dissipation e(r) is very constant for the interval of 
time to < t < tf in the database, varying by < 1% from 
its space-time mean value e = 0.0919 over that interval. 



(a) 



As in our study of Richardson diffusion, we solved (27) 
using the Euler-Maruyma scheme with dt = 10~ 3 and 
solved the equation (44 1 for J~ with the Euler method. 



Velocity-gradients are calculated by a 4th-order finite- 
difference scheme with 4th-order Lagrange interpolation 
in space. We checked convergence in dt for several x 
values both by taking smaller dt and by comparison with 
the 1.5th-order of Platen for ( p7| and with a consistent 
equation for the matrix JT. We then approximated 



(B 2 ( X ,t f )) T = 
2 



]T itr (j(a n , tf ,T)J T (a m ,t f ,T)), 

(46) 



<m<N 



by a sum over the N(N— 1)/2 = 523, 776 number of pairs 
of samples. N ote th at Hoeffding's law of large numbers 
for ?7-statistic d 81 * 82 l implies that this pair-average con- 
verges for N — > oo to the double-average in (45) over 



two independent realizations of the white-noise. Wc then 
furthermore averaged in space over S — 1600 points x s , 
s = 1, 2, ... S, with 25 points chosen randomly from each 
of 64 = 4 3 subcubes of the whole domain. We obtain 



1 S 

(B 2 (i)) = ^5>? 2 (x SJ t/)> 7 

S — 1 



(47) 



as our final approximation to the magnetic energy. More 
space-averaging was required for the kinematic dynamo 
than for Richardson diffusion because of intermittency of 
the velocity-gradients involved in line-stretching^!. 



2. Results 

We now present our results for Pr m = 1. We first 
demonstrate convergence of our algorithm in S and N. 
There are several ways to estimate the errors associated 
with the averaging over space and random samples. One 
approach is to consider the approximation to the Cauchy- 
Green matrix F obtained by omitting the ~ factor and 



the trace in eq.(46). As we have discussed above, the ex- 



act Cauchy-Green matrix should be a positive-definite, 
symmetric matrix. Thus, if we form the symmetric and 
anti-symmetric parts 



1 



(F + F T ), F, 



1 



F T ), 




time t [resistive units] 
(b) 




time t [resistive units] 



FIG. 5. Error estimation, (a) The ratio p(t) of the norms of 
the anti-symmetric and symmetric parts of the approximate 
Cauchy-Green tensor, (b) The eigenvalues <j>%{t), i = 1,2,3 of 
the approximate Cauchy-Green tensor, with red for largest, 
green for middle, and blue for smallest. 



one measure of the relative error in our calculation is the 
ratio of matrix norms 



P(t) 



(FA(tf,T))\\ 

<Fs(*/,t)>|| 



Furthermore, if the small-scale turbulence is statistically 
isotropic (as is known for the database employed), then 
the space-average Cauchy-Green matrix should satisfy 



^tr((F))%. 



In particular, each of the three eigenvalues cf>i(t), i 
1,2,3 of (F(t/,T)) should be equal to (B 2 (t)}. 
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In Fig. 5 (a) and (b) we plot our results for p(t) 
and 4>i(t), i = 1,2,3, respectively. Note that time t in 
these plots and in all those following have been non- 
dimcnsionalized by the resistive time t v = y/\/e = 4.49 x 
10~ 2 (which is also the viscous time since Pr m = 1). The 
results for p(t) show that the relative error in our calcula- 
tion is less than a few percent up until about 15 resistive 
times. This is confirmed by the plot of the three eigenval- 
ues 4>i{t), i = 1, 2, 3, in panel (b), which are in quite close 
agreement until that time. The eigenvalues also remain 
all positive until after 18 resistive times, consistent with 
positivity of (B 2 (t)). These results show, incidentally, 
that there was no need for us to take the initial magnetic 
field Bo to be random and statistically isotropic. The 
same dynamo growth is observed for any deterministic 
uniform field pointing in any direction. 

The error in our approximation for (B 2 (t)} can also 
be estimated in another way. Since the number of pairs 



in (46) is quite large, one can guess that the dominant 



error arises from the average over S space points in (471 



A central limit theorem argument then suggests that the 
error is approximated by 



5{B\t)) c 
with the spatial variance 



Vax{(^(t / )) T )} 
S 



(48) 




FIG. 6. Mean magnetic energy. Plotted in blue is the mean 
magnetic energy (solid line) calculated from (46 1, (47 1, along 



with plus and minus the error (dotted line) estimated from 
( 48 1 , ( 49 1 . The straight line in red shows the least-square lin- 
ear fit over the time interval t = 6 to t = 16. 



Var^i? 2 ^)).)}- 



1 



s 



1 ^ 

s=l 



(B 2 (x S) t f )) T -(B 2 (t)) 



This latter quantity has some independent physical inter- 
est, because it quantifies the spatial intermittency of the 



dynamo effect. In Fig. 6 we plot our approximation (47) 
for (B 2 (t)) along with plus-or-minus the error estimate 
5(B 2 (t)) above. These are consistent with the estimates 
in the previous Fig. 5. 

The plot of magnetic energy from kinematic dynamo 
effect in Fig. 6 is our central result. The expected ex- 
ponential growth ~ e 7 * of magnetic energy is clearly ob- 
served after about 6 resistive times. A linear fit over the 
range of times 6 to 16 is also plotted in Fig. 6, yield- 
ing an estimated growth rate of jt n = 0.158. At ear- 
lier times the growth rate is significantly larger. For 
example, a linear fit over the range of 1 to 4 resistive 
times yields an estimate jt„ = 0.344. This is closer to 
the magnitude of the typical viscous strain rate eigen- 
value y / (S' 2 )/3, which in units of the viscous/resistive 
rate ySjv is equal to l/\/6 = 0.408. The physical in- 
terpretation of these results is clear. According to the 
stochastic Lundquist formula ( |30[ ), field lines that are 
carried into a point x along individual Lagrangian tra- 
jectories are stretched at the viscous strain rate. How- 
ever, the resistive average over the ensemble of stochastic 
trajectories leads to cancellation and suppression of the 
growth rate. As time advances, the spatial region sam- 
pled by the wandering trajectories increases in size and 
the suppression effect increases. Indeed, the asymptotic 




time t [resistive units] 



FIG. 7. Local slope of backward dispersion. Plotted in blue is 
the local slope of the backward dispersion {r 2 (t)) y as defined 
in ( 50 1 . Plotted in red is the horizontal line-segment from 
t = 6 to t = 25 for the Richardson value a = 3. 



exponential growth range begins precisely at the onset 
of turbulent Richardson diffusion of the trajectories. To 
demonstrate this, we consider the (backward) dispersion 
(r 2 (t)) of the stochastic Lagrangian trajectories which 
determines the typical linear dimension L(t) — \/ (r 2 (t)) 
of the region from which field-lines arrive. The local-in- 
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time scaling exponent 



(a) 



a(t) 



dlog((r 2 (t))) 
dlog(t) 



(50) 



will equal 1 at early times and transition to 3 when 
Richardson diffusion sets in. The numerical results for 
this quantity are plotted in Fig. 7, along with a hori- 
zontal line (red) at level a = 3 beginning at t = 6. It 
is clear that the asymptotic exponential growth range of 
magnetic energy and the Richardson diffusion of trajec- 
tories start at the same time. The inertial-range prop- 
erties of turbulent 2-particlc diffusion are thus critical in 
determining the ultimate growth rate of the Pr m = 1 
turbulent kinematic dynamo. 

Although initial field lines arrive to a point at time t 
from a large region of size L(t), not all of these field lines 
contribute equally to the growth of magnetic energy. To 
quantify the contribution of line- vectors at initial separa- 
tion r, we use the line-vector correlation function which 
was proposed as a "dynamo order-parameter"^ 

1l k t(r,t) = Fli{r,t \0,t), 

where F is the homogeneous 2-particle Green's function 



from (42F^. TZ k e(T,t) represents the scalar correlation 
at time t between material line- vectors £(t), £' (t) which 
started as unit vectors , eg at positions displaced by r 
at the initial time and which arrive at the same final 



point. Setting r = in eq.(41 1 gives, in general 



(B 2 (t)) 



d 3 r 



(B*(p)B e (0))K ke (p,t). 



This formula separates the effect of the initial correla- 
tions of the magnetic field and the effect of the turbulent 
advection and stretching. For the case of isotropic and 
non-helical velocity statistics, we may decompose the ten- 
sor 1Z into contributions from line-vectors initially longi- 
tudinal and transverse to the separation vector r : 

fckt(*, t) = R L (r, t)f k r e + R N (r, t)(S ke - f k f e ), 

where f is the unit vector in the direction of r. For our 
particular choice of initial magnetic field, 

1 f°° 

(B 2 (t)) = - / Avr 2 dr [R L (r, t) + 2R N (r, t)]. (51) 
Jo 

This line-correlation can be calculated numerically by 
the same procedure that we have used to obtain the mag- 
netic energy itself in eqs.(46) and (47). For example, the 
longitudinal line-correlation can be approximated by 

R L (r,t) = 



ti N ( N - 1) : 



, «*(|a?(r)-ar(T)|-r) 
4irr 2 




r resistive units 



(b) 




[resistive units] 



FIG. 8. Line-vector correlations: (a) longitudinal, (b) trans- 
verse. The line- vector correlations Rl(t, t) and fly(r , t), re- 
spectively, as calculated numerically from formula ( |52[ ) at 
time t — 11.12, are plotted in blue with solid lines and the 
correlations plus and minus their estimated errors with dot- 
ted lines. Also plotted with red lines are analytical resulti^ 
for these correlations in the Kazantsev-Kraichnan model at 
Pvm = 0, normalized as described in the text. 



taking f as a row vector. The corresponding transverse 
correlation Riq{r^t) is obtained by replacing f in (52) 
with two orthogonal unit vectors ej, i = 1,2 which span 
the subspace orthogonal to f and by then summing over 
these two contributions. In practice these continuous dis- 
tributions must be sampled in discrete bins. We took 200 



(52) 



bins of size Ar = l n /2, or one-half of the resistive length- 
scale. To capture the contributions from r > I00£ v but 
to avoid large fluctuations in the results, we added three 
extra large bins of size 100, 200, and 400 £ v , centered at 
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KK model with spatial Holder exponent h = 2/3 in the 



20 40 60 80 100 120 140 160 180 200 

r [resistive units! 



FIG. 9. Contributions to magnetic energy. Plotted in 
blue is the function 4nr 2 Ri,{r,t) and in cyan the function 
4-irr 2 R]s/(r,t), both at time t = 11.12. These represent the 
contributions to magnetic energy in (51 1 from pairs of line- 



vectors at initial separations r and initially parallel and per- 
pendicular, respectively, to the separation vector r. 



150,300 and 600 £ v , thus covering the whole range of r 
in the database. 

We present in Fig. 8 our numerical results for Rl{t, t) 
and Rjy(r,t) at a time t = 11.12 in resistive time-units. 
As should be clear from Figs. 6 and 7, this time lies well 
within the range both of exponential growth of energy 
and of Richardson diffusion of particle pairs. The most 
interesting feature of these correlations is their consid- 
erable diffuseness in r. This is shown even more clearly 
in Fig. 9, which plots the integrands 47rr 2 i?i(r, t) and 
4nr 2 R N (r,t) in eq.|5l]) for (B 2 (t)). In order to get 50% 
of the magnetic energy one must integrate in that formula 
out to r = 15. 71£rj- Likewise, to get 75% of the energy 
one must integrate out to r = 39.15^ and to get 90% 
one must integrate all the way out to r = 66.59^. Thus, 
line-vectors separated initially by many resistive lengths 
are brought together by turbulent advection to produce 
the dynamo growth. Another very interesting feature is 
the large negative contribution of the initially longitudi- 
nal line-vectors, seen in Fig. 8(a) and even more clearly in 
Fig. 9. A very similar negative contribution was found 15 
for the Kazantsev-Kraichnan dynamo model at Pr m = 
and Re m = oo. It was suggested there that negative val- 
ues of Rl (r, t) for sufficiently large r are due to an effect 
of bending and looping of field lines. 

It is interesting to make a more detailed compari- 
son of our numerical results for hydrodynamic turbu- 
lence at Pr m = 1 and of the analytical results in the 
Kazantsev-Kraichnan (KK) dynamo model for Pr m = 0. 
The comparison is not completely straightforward, so a 
few words of explanation are required. We consider the 



velocity correlation (15 1. With this choice of h, Richard- 
son's i 3 -law holds for particle dispersion and it is gener- 
ally the appropriate choice for comparisons with physi- 
cal fluid turbulence. In this case, the coefficient D\ in 
(15) has the same dimensions as e 1 / 3 , with e the mean 



energy dissipatio n per unit mass. Our previous results 
for the KK modefHUS were obtained with lengths non- 
dimensionalized by (A/-D1) 3 / 4 and times by (A/JD 3 ) 1 / 2 . 
However, D\ is not numerically equal to e 1 / 3 , but only 
dimensionally the same. In order to relate them quan- 
titatively, we note that the KK particle diffusion equa- 
tion ( 16 ) reduces in the isotropic sector to Richardson's 
original equation ^ with coefficient of cddy-diffusivity 
K a = 2D 1 . In that case, one recovers the t 3 -law ([6j with 
a Richardson-Obukhov constant go if one makes the def- 
inition for the KK model 



64 o 
143 £>?. 

8lg 1 

1/4 



We thus obtain £ v = (A 3 /e 



8 = -A= i- 35 - 

2V2 V 143 



,1/4 



(53) 

/3(A/L>i) 3/4 with 
For the backward 2-particle diffu- 
sion in our kinematic dynamo study we have found that 
go = 1-57 (somewhat larger than the value 1.35 reported 
in section II. C for an independent experiment over a dif- 
ferent time range) and thus 8 = 0.343. For comparison 
with the present results, therefore, the results of Eyinlfl^, 
Fig. 10 must have £-axis scaled by 1/(3 and j/-axis scaled 
by /3 3 . Furthermore, the quantities Gx(r), Gj^{r) previ- 
ously calculated^ are dynamo growth eigenmodes which 
dominate the behavior of Rl{t, f), Rn{?i t) a t long times. 
Since there is an arbitrariness in the normalization of the 
eigenmodes, we have additional freedom in the vertical 
scale. This is fixed by imposing on the KK eigenmodes 
the normalization condition (51 1 at t = 11.12 resistive 
time units. 

We have plotted the results for the line-correlations 
from the KK model in Fig. 8 (red lines), together with 
our numerical results for hydrodynamic turbulence (blue 
lines). Clearly, the two sets of results are qualitatively 
very similar. In both cases, the transverse correlation 
iZjv(r, t) is everywhere positive, sharply peaked at small 
r, but with a slow decay at large r. The longitudinal 
correlations Rl (r, t) of the two sets share all these same 
features except for their sign, with both exhibiting a long 
negative tail at large r. The major difference between the 
results for hydrodynamic turbulence at Pr m = 1 and the 
KK model at Pr m = is the distinctly slower rate of 
decay of correlations at large r for the latter. 

The precise decay rate at large r is known for the line- 
correlations in the KK model. It was shown 15 that, up to 
power-law prefactors, these exhibit a very slow stretched- 
exponential decay 



G L (r), G N (r) ~ exp - 



2L>i 



/3 



r > I,,, (54) 



where 7 is the dynamo growth rate. A remarkable fea- 
ture of this asymptotic formula is that it depends upon 
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(a) 



stretched exponential. If we use the relation ( 53 1 to re- 



-turbulence : 
- KK model 
asymptotics - 




1 



r 1 / 3 [resistive units] 
(b) 




FIG. 10. Line-correlations for large separations: (a) longitu- 
dinal, (b) transverse. Plotted in blue are the numerical results 
from (52 1 at time t — 11.12 and in red the analytical results^ 
from the Kazantsev-Kraichnan model at Pr m — 0. Straight 
lines in the log-linear plot versus r 1//3 correspond to the large- 
r asymptotics (541. Shown in is the prediction of that 

formula with the growth rate 7 = 0.158 determined in Fig. 6. 



resistivity rj (or magnetic diffusivity A = r/c/iir) only 
through the growth rate 7. To check for a similar de- 
cay in hydrodynamic turbulence we present in Fig. 10 
(a),(b) a log-linear plot of the line-correlations |i?z,(r,i)|, 
i?./v(r, t) versus r 1 / 3 , both for our numerical calculation 
with the hydro turbulence database (blue) and for the 
KK model (red). The straight lines in these plots ver- 
ify the stretched exponential decay with power r 1 / 3 . In 
fact, the decay law (54) which was derived in the KK 
model for Pr m = holds very well in hydrodynamic 
turbulence at Pr m = 1, including the coefficient in the 



place Di with e and then substitute the dynamo growth 
rate jt^ = 0.158, we predict from (54) a line in a log- 
linear plot with slope —3.507 (in resistive units). A line 
with this slope is plotted (in green) in Fig. 10 (a),(b) and 
can be seen to match our numerical results quite well. 
The reason for the slower decay of the stretched expo- 
nential in the KK model at Pr m — is thus entirely 
attributable to a smaller dynamo growth rate than what 
we find in hydro turbu lence for Pr m = 1. In particu- 
lar, using the resulP^ that ^{X/Df) 1 / 2 = 0.193 for KK 
at Pr m = and using again the relation (53), we get 



that jt n = 0.0217. This is almost an order of magnitude 
smaller than the result 7< r) = 0.158 that we found for hy- 
dro turbulence at Pr rn = 1 and motivates some further 
discussion below. 



3. Discussion 

Let us make a quantitative comp arison of ou r growth 
rate with those found in other studies^ 21 * 72 * 73 * 85 -] of small- 
scale dynamo in incompressible, non-helical turbulence 
at Pr m = 1. We present both the Reynolds numbers 
and the dynamo growth rates found in these studies. We 
give the "box-size" Reynolds number Re = v! /vk$, where 
u' = u rms /\/3 is the rms value of a single velocity compo- 
nent and kg is the smallest wavenumber in the simulation 
study. This is not as dynamically significant a Reynolds 
number as is Ref = u rms /vkf, based on the forcing 
wavenumber kf, but it is the easiest to calculate from the 
published data. We give the growth rate in the dimen- 
sionless form r yt v where, as above, t v — {v/e) 1 / 2 is the 
viscous time (and also the resistive time for Pr m = 1). 
The data from the paper of Schekochihin et al. (2004 J 3 -' is 
taken from their Table 2 (where llty = vToT rms there). 
The data of Haugen et al. (2004JFU is taken from their 
Fig. 3, which plots their results for a = (7 / '2) / \u rms k f) 
versus Ref. Note, however, that 



jt u = 2a 



3i?e T 



5 Re 



f 



where Rex = u'£t/v is the Taylor-scale Reynolds num- 
ber defined there^H with It = V^u rms /u) rms . We can 
also infer from their Table iP^ that Rex = \J 60i?e / so 
that, putting all these relations together, 



12a 



Finally, Re = (1.5) Rej since their forcing wavenumber^ 
is kf — 1.5. The paper of Vincenzi (2002)P3 presents a 
numerical study of the KK model. Fig. 6 in that work 
plots the results for ^(X/Df) 1 ' 2 at Re — 00 and a broad 
range of Prandtl numbers. Note that the "viscosity" for 
the KK model is definecP^ by = Diit^ where l v is 
the short-distance cut-off length of the inertial scaling 
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range and the corresponding "Prandtl number" is de- 
fined by Pr mie = i/*/A. It is not obvious how to best 
compare this "Prandtl number" with that for viscous 
hydrodynamic turbulence. In any case, the result for 
Pr m * — 1 is 7(A/Z?f) 1 / 2 = 0.350, which we transform 
to 7t„ = 0.0383 as described earlier. Finally, the papers 
of Cho et al! 72 * 73 * present results for a large number of 
runs with both normal and hyperviscosity. The authors 
have kindly provided me with data from their highest 
Reynolds-number normal-viscosity simulation, denoted 
RUN 512P - -B 10~ 3 (for 512 3 resolution with physical 
viscosity and initial uniform seed field of strength 10 -3 .) 
This data is reported below (J. Cho, private communica- 
tion). In one important respect this simulation is quite 
distinct from the others reported here, because the ini- 
tial velocity field is not a fully developed turbulent field 
but is instead supported in the very low-wavenumber in- 
terval 2 < k < 4. Turbulence quickly develops with the 
volume- averaged kinetic energy dissipation e(t) increas- 
ing by a factor of 42 over their time interval of kinematic 
dynamo. Although this dynamo simulation is not in a 
statistically steady turbulent regime, it does give some 
useful perspective and we thus include it here. All of 
these data are gathered into Table I below. 



TABLE I. Reynolds numbers and dynamo growth rates from 
several numerical studies with Pr m = 1. Shown are the box- 
size Reynolds number Re — u' /i/ko, where ko is the minimum 
wavenumber, and the growth rate 7 non-dimensionalized by 
the viscous time t u = (u/e) 1 ^ 2 . 



Reference 


Re 


7t„ 


Schekochihin et al. (2004) 


110 


0.0138 




210 


0.0265 




450 


0.0329 


Haugen et al. (2004) 


227 


0.0338 




450 


0.0379 




1050 


0.0385 


Vinccnzi (2002) [KK model] 


00 


0.0383 


Cho et al. (2009) 


294 


0.550 


present study 


6380 


0.158 



It is widely expected that r yt u approaches a univer- 
sal value as Re — > 00, but the data presented seem not 
to confirm this. The study of Haugen et alPS indeed re- 
ported seeing such a limiting behavior, as can be observed 
from their data in Tabic I. Their results for the growth 
rate are also consistent with those of Schekochihin ct aP 
and remarkably close to those of VincenzPSl f or the KK 
model at Re = 00. The latter must be regarded as a 
coincidence, however, due to the ambiguity in the defini- 
tion of the Prandtl number for that model. For example, 
a better definition might be Pr m = e 1 / 3 ^!/ 3 / A, which, 
using (53), gives Pr m = 4.16Pr m *. Thus, Pr m = 1 cor- 
responds to Pr m * = 0.24. Whatever the "best" definition 
of the Prandtl number for the KK model might be, the 
value of jt,, for Pr m = 1 will be somewhat smaller than 



that reported above for Pr m * — 1, because the latter is 
near the maximum of 7(A/D 3 ) 1 / 2 (VincenzP^, Fig. 6). 
Our own result for r yt v is about 4 times larger than the 
value of Haugen et al. for their highest Reynolds number 
simulation. The dimensionless growth rate increases with 
Re and our Reynolds number is about 6 times larger than 
theirs, so this might account for some of the discrepancy. 
On the other hand, our resolution of the viscous range is 
also relatively poor and this may degrade the numerical 
accuracy of our result for the growth rate. To further 
complicate the picture, the value of jt u found by Cho et 
al. is about 3.5 times larger than ours. Note that the 
result reported in Table I for Cho's data used the value 
e from a time-average of e(t) over the kinematic interval 
of exponential growth in order to define f „ . If we instead 
used the value e(i*) at the end of the kinematic interval 
to define i„, then we would obtain a somewhat smaller 
value jt u = 0.312 for Cho et al. but still twice ours. 

The results in the Table I appear contradictory on the 
face of it, but we believe there is a simple resolution. The 
key is the time range over which the exponential growth 
is observed. Our growth rate was calculated for the time 
interval 6-16^ with the latter time only about 1/3 T Ul 
where T u — L u /u' is the large-scale eddy-turnover time. 
Haugen et al. do not publish growth curves in their pa- 
per, but Schekochihin et al. present in their Fig. 
a log-linear plot of the magnetic energy for their kine- 
matic runs, with fitting ranges for exponential growth 
indicated. The cited growth rate for their 450 Reynolds- 
number simulation (Run A) was obtained over a time- 
interval of 66. 6-505. 9t„. This corresponds in their sim- 
ulation to about 5.86-44.51^. (Here we have used the 
well-known result that T u = 0A(u') 2 /e and the data in 
their Table P to estimate T u = 0.267 for Run A.) At 
earlier times in their simulation, the growth of magnetic 
energy is much faster. For example, if we use the data in 
their Fig. 21 for times < 5.86T U to estimate the growth 
rate, we obtain a value yt y — 0.0691, about half of ours. 

We therefore conjecture that there are two distinct 
kinematic regimes with exponential growth of magnetic 
energy at different rates, one for times t <C T u and an- 
other for times ( > T„. This makes perfect sense from 
the Lagrangian point of view developed in this paper. 
For times t <C T u , the (backward) particle dispersion 
(r 2 (t)) is growing as ~ et 3 (see Fig. 7). This means 
that the magnetic field strengths are obtained by av- 
eraging initial field vectors arriving from inertial-range 
separations growing as a power-law. However, for times 
t^>T u , the mean-square separation (r 2 (t)) saturates to a 
value of order ~ L 2 B , where Lb is the size of the periodic 
box. In this time regime, magnetic field strengths are 
obtained by averaging seed field vectors arriving to each 
point from the entire flow domain. It stands to reason 
that the growth rate would be reduced in this strongly 
mixed regime. Which of the two kinematic regimes is 
most relevant in practice depends upon t he stre ngth 
of the seed magnetic field B . It is observe d 3 * 21 * 72 ! that 
the kinematic interval of exponential growth ends once 
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the magnetic field energy (B 2 (t)) grows to equipartition 
with the viscous-range kinetic energy u 2 = {ev) 1 / 2 . If 
]n(u„/B a ) < jT u = Oi(Re) 1 / 2 ), then this occurs at 
times much less than T u , the regime that we have stud- 
ied. If the seed field is extremely small or the Reynolds 
number not so large, however, then viscous equipartition 
will only occur for times 3> T u . This is the regime studied 
in Shekochihin et alP and Haugen et alPS In support of 
this conclusion, we note that the RUN 512P — PolCU 3 of 
Cho et al., with a not too small seed field, reached satu- 
ration with viscous kinetic energy at a time around 12^ 
or 1.2T U . The growth rate cited in Table I was obtained 
from a fit over the range 1.72-11. 7t„ or 0.173-1. 18T U , 
comparable to the time range that we study, and with a 
similarly large growth rate. 

Clearly more study of this matter would be desir- 
able. Unfortunately, the turbulence database that we 
employ only stores the velocity field for about one T u , or 
44.8i„, so we cannot study times ^> T u . We do believe, 
however, that our results have attained the asymptotic 
regime t v <C t <C T u . In support of this claim, we note 
that we are observing already at time 11.12i„ the ex- 
pected inertial-range behaviors, such as Richardson dif- 
fusion and the stretched-exponential correlation decay in 
eq. ( 54 ) (which was derived in the KK model for L u = oo) . 

The most important general conclusion from our re- 
sults presented above is that the inertial-range phe- 
nomenon of Richardson diffusion plays a critical role in 
the small-scale fluctuation dynamo of hydrodynamic tur- 
bulence at Pr m — 1. Because field lines are only "frozen- 
in" stochastically, an infinite number of lines enters each 
point from a very large region of size L(t) ~ (et 3 ) 1 ^ 2 
at time t. This mixing of field lines from far away op- 
poses the growth by stretching of the individual lines and 
suppresses the dynamo growth rate. Because of nearly 
complete cancellation, the lines arriving from separations 
of order ~ L(t) contribute a vanishingly small amount 
to the magnetic energy. Nevertheless, the asymptotic 
formula (54) shows that field- lines separated at inertial- 



range distances r>l, give a non-negligible contribution. 
That formula implies also — very remarkably — that the 
dynamo growth rate can be inferred directly from the 
stretched-exponential decay of the magnetic line-vector 
correlations in the inertial range. 

A second important conclusion is that there are very 
fundamental similarities between the mechanisms of the 
small-scale fluctuation dynamo for Prandtl numbers 
Pi"m — 1 & n d Pr m — in high Reynolds-number tur- 
bulence. The physics that we have discussed above for 
hydrodynamic turbulence at Pr m = 1 an d larg e Re m 
agrees very closely with what was establishe d 14 ! 15 ! for the 
Kazantsev-Kraichnan model at Pr m = and Re m = oo. 
Richardson diffusion and stochasticity of flux-freezing 
play a critical role in both. We thus disagree with 
Shekochihin et alP who drew the conclusion that "The 
small-scale dynamo at Pr m ~ 1 belongs to the same class 
as the large-Pr m dynamo." There is a flaw in the ar- 
gument made by those authors that inertial-range time- 



scales are too long to affect the operation of the dynamo, 
which proceeds at a faster viscous rate (see quote in the 
opening paragraph of Section IV. B). Because the expo- 
nential dynamo growth is a long-time phenomenon, there 
is sufficient time for advection by inertial-range eddies to 
affect and modify its rate. Indeed, we have seen in the 
present study at Pr m — 1 that the inertial-range phe- 
nomenon of Richardson diffusion commences at times 
t > &t v and not at times orders of magnitude greater 
than the viscous time. 

In our view, the Pr m — 1 small-scale dynamo — and, 
more generally, the finite Pr m dynamo — is a transitional 
case which shares some of the attributes of both of the 
extreme limits Pr m = and Pr m = oo. A corollary 
of this view is that there should be an important influ- 
ence of the inertial range even for small-scale dynamos 
with Pr m ^> 1, if also Be> 1. This is the situation for 
many low-density, high-temperature astrophysical plas- 
mas with large magnetic Prandtl numbers but with, also, 
a substantial inertial range. As we have discussed in sec- 
tion II. B, the exponential separation of particles typical 
of the large- Pr m "Batchelor regime" is an intermediate 
asymptotics for a range of times i„ -C t <C \n(Pr m )t v 
in high Reynolds-number turbulence. After a time t > 
t„ln(Pr m ), magnetic line-elements in such a turbulent 
flow will begin to experience Richardson diffusion and a 
concomitant decrease of the dynamo growth rate. Be- 
cause the dependence upon Pr m is logarithmic, the time 
to enter this inertial-range influenced regime is only a 
relatively small multiple of t v even at very large Pr m . 
(For example, with a value Pr m = 10 14 typical of the 
warm interstellar medium, ln(Pr m ) = 32.) If the initial 
magnetic seed field is weak enough, then its backreaction 
on the turbulence up to that time may be neglected and 
the turbulent kinematic dynamo process essentially as we 
have discussed in this section of our paper will proceed, 
even at very large Prandtl numbers. 



V. NONLINEAR MHD TURBULENCE, DYNAMO AND 
RECONNECTION 

Although we have focused in our discussion of dynamo 
effect on the kinematic stage, almost all of the results 
of this paper extend — with appropriate modifications — 
to fully nonlinear MHD turbulence with backreaction on 
the flow from the Lorentz force. A full treatment is not 
possible here, but we shall try to stress what is general 
in our previous presentation, sketch any necessary modi- 
fications, discuss relevant references, and point out some 
important directions for further work. 

The phenomenon of 2-particle Richardson diffusion 
and "spontaneous stochasticity" will doubtless exist in 
MHD turbulence. Their cause is the roughness of the ad- 
verting velocity field and all theories, simulations and ob- 
servations of MHD turbulence (whatever their other dif- 
ferences) agree that both the velocity and magnetic fields 
are indeed rough in the inertial-range of such flows. The 
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quantitative growth law of particle dispersion will depend 
upon the precise characteristics of MHD turbulence, such 
as the spectral slopes, degree of anisotropy, etc. These 
depend upon the ultimate theory of MHD turbulence, 
which is an open problem. We base our discussion be- 
low on the Gold reich-Sridhar (GS) theory of strong MHD 
turbulenc e! 86 * 87 !, but alternative theoriegS^D would lead 
to similar results. 

For simplicity, we assume incompressible and sub- or 
trans- Alfvenic MHD turbulence, with rms velocity fluc- 
tuations u' < v A , where v A = B / ' y/Airp is the Alfven ve- 
locity based on the external magnetic field strength Bq. 
We also assume that large-scale anisotropy is such that 
"critical balance" holds throughout the inertial range and 
the turbulence is strong. Finally, we assume that mean 
cross-helicity is zero and there is an equal flux of upward 
and downward propagating Alfven waves. Under these 
conditions, the GS theory predicts that velocity incre- 
ments for separations I scale as 

5u{i) ~ (el!.) 1 / 3 ~ {el\\/v A ) l /\ 

where e is energy dissipation per mass and £±, l\\ are the 
separations perpendicular and parallel to the magnetic 
field, respectively. We shall also assume that perpendic- 
ular (shear- Alfven) and parallel (pseudo- Alfven) compo- 
nents of the velocity increments scale in the same way, 
or 5uj_(£) ~ 5u\\(£). In that case, we can give a simple 
theory of particle-pair dispersion like that presented for 
hydro turbulence in section II. A. If we let r±(t) and ry (f) 
be the particle separations perpendicular and parallel to 
the field, respectively, then 



dt 



(assuming that rn = 0) implies that 



et 6 



(55) 



and in the same manner 
d 



dt 



r\\ ~ 5u\\{r) ~ {er\\/v A ) 1/2 



(assuming that r± = 0) implies that 
r\{t) ~ {e/v A ft\ 



(56) 



The above scaling laws will hold for intermediate times 
where t > (r^^/e) 1 / 3 , (v A r\\ (0)/e) 1/2 but t < T u = 
L u /vf, the large-eddy turnover time. The most impor- 
tant feature observed in both cases is that initial sep- 
arations are "forgotten", the physical basis of the phe- 
nomenon of spontaneous stochasticity. A full theory 
of 2-particle dispersion in MHD turbulence will obvi- 
ously be quite intricate and will depend upon the regime 
of turbulence considered and the phenomenological as- 
sumptions e mp loyed. Within theories of weak MHD 
turbulencc422H2il ; it is should be possible to give an ana- 
lytical treatment using well-established methods^. How- 
ever, whatever final form the theory of MHD turbulence 



may take, spontaneous stochasticity seems to be a likely 
consequence. 

There have been a few numerical studies of 2- 
particl e disp ersion in MHD turbulence, both wittP^ and 
withoulpZESl external magnetic field. In the presence of 
an external magnetic field, Busse and Miiller find that 
r±(t) grows faster than r\\(i) (see their Fig. J^ID), in quali- 
tative agreement with our formulas (551 and (56l 9 ^. The 



situation is not entirely clear, however. Following Cho 
and VishniacP^, it is generally believed that similar align- 
ments and anisotropics will hold at small scales in MHD 
turbulence without an external magnetic field, just as for 
the external field case, if the alignments are taken with 
respect to a local magnetic field. Paradoxically, however, 
Busse et al. find in MHD turbulence without external 
field that r±(t) <C rn(t) when these quantities are defined 
with respect to the local field. The total displacement 
vector r(t) thus becomes preferentially aligned with the 
local magnetic field. It is not obvious how to explain this 
observation, but it possibly has something to do with 
the dynamical alignment of velocit y and magnetic-field 
increments at small scale d 91 * 100 * 101 *. Clearly, more study 
of these issues is required. 

Our results in section III. A on stochastic flux freez- 
ing in resistive magnetohydrodynamics were derived in 
a fully nonlinear setting. No kinematic assumption was 
made there. It has elsewhere been showrp2l that all ef- 
fects of the Lorentz force on fluid motion are described 
by a second stochastic conservation law which gener- 
alizes the Kelvin circulation theorem, at least for in- 
compressible fluids with Pr m = 1. This result extends 
to resistive MHD th e "gene ralized Kelvin theorem" de- 
rived for ideal MHD^JiMI. As we shall show in a fu- 
ture paper, this additional stochastic Kelvin theorem also 
holds in compressible plasma fluids with Pr m = 1 if 
they are barotropic (pressure depending only on density 
and not on temperature) and in non-isothermal fluids 
if the thermal Prandtl number is also unity. The ex- 
istence of two stochastically "frozen-in" fields provides 
strong constraints which deserve to be further explored. 
We remark finally, and most importantly, that our dis- 
cussion in section III.B on the high-Reynolds limit of 
plasma turbulence made no kinematic assumption. The 
prediction that stochastic flux-conservation holds in that 
limit depends only upon the phenomenon of spontaneous 
stochasticity, which we have argued should inevitably oc- 
cur in high Reynolds-number plasma turbulence. 

The Lagrangian formalism of magnetic dynamo pre- 
sented in section IV. A is valid in nonlinear MHD turbu- 
lence, if averages over velocity ensembles are conditional 
on a given initial magnetic field. An interesting ques- 
tion to pursue, even in the kinematic stage, is the role 
of compressibility on the small-scale, turbulent dynamo. 
The behavior of Lagrangian particles in strongly com- 
pressible flows are dramatically different from those in 
incompressible flows and this should have a signficant 
impa ct on d ynamo action, e.g. possible Prandtl-number 
effects^nEH. Our La grangian numerical methods in sec- 
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tion IV. B are specific to the kinematic dynamo. However, 
it would be possible, and interesting, to employ them to- 
gether with simple phenomenological model s of nonlinear 
saturation that have been propose d 104 * 105 !. This would 
obviate the need for additional closure approximations 
and would help to understand some of the physics of the 
saturation process. 

Perhaps the most important implications of the present 
work are for the problem of turbulent magnetic reconnec- 
tion. Our results and arguments show compellingly, we 
believe, that the constraint of flux-freezing in a turbu- 
lent plasma at high-conductivity must be quite differ- 
ent than is generally understood. The naive estimate 
of flux-line slippage due to resistivity, (r 2 (t)) ~ At, is 
incorrect, by many orders of magnitude. The correct es- 
timate will depend upon the ultimate theory of MHD 



turb ulen ce, but it must have the general form of ( 55 ) 



and (56) above. The quantity \J (r^_(t)) can be inter 
preteaas the lateral distance that magnetic field-lines 
diffuse through a turbulent plasma in time t, and ( 55 ) 



should be particularly useful in estimating reconnection 
rates for astrophysical phenomena. Note that both (55) 



and ( 56 1 are completely independent of the resistivity 
(and of any other microscopic plasma mechanism of line- 
slippage). These results give s upport, therefore, to all 
theories^SHlMI anc j obser vat ions^QJlLEl f f as t magnetic 
reconnection in MHD turbulence. In fact, our results 
show that fast reconnection is necessary for (implied by) 
all standard theories of MHD turbulence^HSU. 

There are particularl y close connections with the 
Lazarian-Vishniac theory^ 7 * 108 ! based on stochastic wan- 
dering of field-lines. We have stressed that Lagrangian 
particle trajectories become intrinsically stochastic in a 
rough velocity field, due to the phenomenon of Richard- 
son 2-particle dispersion. The "stochastic wandering" in- 
voked by Lazarian-Vishniac is an analogue of Richardson 
diffusion for the field lines themselves in a rough magnetic 
field. Because of this effect, there is not just one field-line 
passing through each point in the limit of zero resistiv- 
ity, but instead an infinite ensemble of random field-lines. 
The con nection between our ideas and those of Lazarian- 
Vishniac 107 108 deserve to be further examined. 



VI. FINAL DISCUSSION 

Two ideas are commonplace in the literature on plasma 
magnetohydrodynamics. One is that flux-freezing must 
hold approximately for Re m 3> 1, an assumption widely 
employed in treatments of magnetic dynamo effect. An- 
other idea frequently advanced in discussion of magnetic 
reconnection is that the flux-freezing constraint is bro- 
ken by rapid diffusion of field lines across thin current 
sheets. There has long been a tension between these two 
ideas, never fully reconciled. Both ideas are sometimes 
invoked in the same setting. For example, it is recog- 
nized that while stretching of "nearly" frozen-in lines 
drives magnetic dynamo action, nevertheless fast mag- 



netic rec onnection i s necessary to relieve tangled field-line 
structur e 1 108 * 115 * 1 ^ . We have argued that both ideas are 
correct, if suitably understood. Magnetic flux through 
individual material loops, advected by the plasma in 
the usual sense, will not be conserved for Re m — > oo, 
because of the development of singular current sheets 
and vortex sheets. However, magnetic flux will never- 
theless be conserved on average for a random ensemble 
of loops, in a novel statistical sense associated to the 
spontaneous stochasticity of Lagrangian flows for high- 
Reynolds-number turbulence. 

Spontaneous stochasticity is a fluid-dynamical phe- 
nomenon due to the explosive separation of particles pro- 
duced by turbulent Richardson diffusion. Because of this 
effect, fluid particles that start infinitesimally close to- 
gether will separate to a finite distance in a fixed amount 
of time, independent of the Reynolds number. In the 
limit Re oo there is an infinite ensemble of Lagrangian 
trajectories for each initial particle position. There is 
already good evidence for this effect from turbulence 
simulations and laboratory experiments, much of which 
is reviewed in Section II. C together with own simula- 
tion results. We expect confirmation from future stud- 
ies at higher Reynolds numbers. The Gottingen Turbu- 
lence Tunnel now in operation should reach Taylor-scale 
Reynolds numbers Rex ~ 10 4 and turbulent Richardson 
diffusion shall be one of the main subjects of investigation 
(E. Bodenschatz, private communication). Richardson 
diffusion and spontaneous stochasticity should also occur 
in high-Reynolds number MHD turbulence and can be 
studied by both simulation and experiment. Recent lab- 
oratory studies of magnetic dynamo and tur bule nt induc- 
tion with low Prandtl-number liquid metalJ^^^ have 
Reynolds numbers high enough and inertial-ranges suf- 
ficiently extensive to support the phenomena. Although 
extremely challenging, it would be quite informative to 
develop particle-tracking techniques for such flows which 
could investigate Lagrangian mechanisms. 

In an earlier work-^H we have considered Alfven's the- 
orem for turbulent plasma flows from a complementary 
perspective. The approach used there was a spatial 
coarse-graining of the MHD equations over a continu- 
ous range of scales £, similar in spirit to a renormaliza- 
tion group analysis. The effective induction equation at 
length-scale I contains a "turbulent EMF" induced by 
subscale plasma motions. As long as energy remains fi- 
nite for high kinetic and magnetic Reynolds numbers, 
the turbulent EMF can be shown to be much larger 
at inertial-range length-scales t than viscous and resis- 
tive terms, or than other possible microscopic dissipa- 
tion terms. In that case, the effective flux-conservation 
equation at inertial-range length-scales I takes the form 



d 
dt 



(p Af(x,i)-dx= <t> £f(x, i)-dx, 



(57) 



where A^ is the coarse-grained magnetic vector-potential 
and x.£ (a, t) is the Lagrangian flow map generated by the 
coarse-grained velocity field up. For a smooth, laminar 
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solution of the ideal MHD equations, eg — ¥ everywhere 
in space as I — > and flux-conservation in the standard 
sense is recovered. However, it was showrP^U that the 



righthand side of (57 1 need not vanish for very singu- 



lar velocity and magnetic fields, in particular, those with 
(coincident) current sheets and vortex sheets. If such se- 
vere singularities appear in MHD turbulence at infinite 
Reynolds numbers, then Alfven's theorem in its conven- 
tional sense can break down. Magnetic flux $t(C, f) = 
$xt(c t) -^(x, t)'dx need no longer be conserved for indi- 
vidual loops C as I — > 0. Here we extend that conclusion 
by arguing that flux-conservation will hold on average for 
a suitable random ensemble of loops. Similarly as in sec- 
tion III.B, we may consider the set of loops a>(C + eC, t) 
at time io obtained by adding random perturbations of 
size e to C and then advecting backward in time with 
a>(-, t) = x7 (•, t). We expect that this ensemble of loops 
in the limits first of vanishing magnetic diffusivity and 
viscosity A, v — > 0, then £ — > and finally e — > will ex- 
actly coincide (at least for incompressible flow) with the 
ensemble a A (C, t) for A, v — > considered in section III. A. 
If so, then flux-conservation will hold in an average sense 
similar to ( 33 1 . The advantage of the present argument 



is that it makes no explicit reference to Spitzer resistivity 
or to any other microscopic plasma mechanism. Stochas- 
tic flux-freezing is fundamentally a phenomenon of the 
nonlinear MHD dynamics. 

Stochastic flux-conservation is expected to be a prop- 
erty of singular solutions of the ideal MHD equations 
that describe turbulent plasmas asymptotically at high 
kinetic and magnetic Reynolds numbers. The existence 
of such singular solutions is an old conjecture of Lars 
Onsager^SHSU Such solutions must be quite different 
from smooth laminar solutions of the ideal MHD equa- 
tions that are familiar from current analysis, however, 
and must possess many "strange" properties. We are 
still learning how to deal with such solutions. We have 
recently shown how to derive the stochastic Kelvin theo- 
rem for the incompressible Navier-Stokes equation using 
a stochastic least-action principle^^l. It is expected that 
there will be a similar stochastic least-action principle 
for Onsager's singular Euler and ideal MHD solutions. 
This is a fundamental motivation to expect stochastic 
flux-conservation at infinite Reynolds numbers. 

There must be many important applications of 
stochastic flux-freezing in plasma physics, astrophysics 
and geodynamo studies. We have here presented one 
concrete application to the finite- Pr m , kinematic, fluctu- 
ation dynamo in non-helical, incompressible fluid turbu- 
lence. Stochastic flux-freezing is critically important to 
the mechanism of small-scale turbulent dynamo, because 
distinct field-lines that are initially separated by inertial- 
range distances arrive to the same point and resistively 
merge to produce the net magnetic field. See Eq.(54). 



Our results and analyses point to essential similarities 
between the finite-Pr m and Pr m = dynamo in their 
Lagrangian mechanisms. Understanding the saturation 
effect of the Lorentz force from a Lagrangian perspective 



is an obvious next step. In addition, there will be in- 
teresting applications of stochastic flux-freezing to many 
other problems, e.g. to the theory of fast turbulent re- 
connection. This will be the subject of future work. 
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Appendix: Path-Integral Formulas 

We here sketch the derivation of the path-integral for- 
mulas ([9]) and (11) in the text. There are discussions 
already available in the literatur e 128 ! 29 ^ , but we stress here 
some connections not found in those works, with rigorous 
stochastic analysis, on the one hand, and with Feynman 
path-integral methods, on the other. 

Our starting point is the SDE ^ , which we discretize 
using the Euler-Maruyama scheme: 



X„ = X n _i + u(X n _ 1 ,t„_ 1 )^ + V2k(W„ - W„_i). 

(A.l) 

Note that the Brownian-motion variables at the discrete 
set of times have the Gaussian density (with Wq = 0) 



V(W 1 ,...,W N ) = 
(const.) exp 



\ n— 1 



W„_i| 



(A.2) 



with respect to DW = Y[ n (£W n . We can obtain the den- 
sity V(X.) from the change of variables formu la "P( X) = 



P(W)/|det (f^)| ■ It is easy to see from (|A.l|) that 



the (3A) x (3N) Jacobian matrix J?^ is block lower 
triangular, with diagonal blocks 4^ 



<9W„ 

ax ■ 



/2/cIj for the 



3x3 identity matrix I. Thus, det ( J^) = const, and 



P(Xi X W )« 

/ 1 N 

\ n=l 



X„_i 



St 



u(X„_i, t n -i 



(A.3) 
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Integrating this density with respect to DX — Y[ n dX n 
and taking the continuous-time limit St — > 0, N — > oo, 
one formally obtains formulas like ^ and (111. 

The mathematically rigorous versions of these 
path-integral for mulas is the classical Girsanov 
transformatioi^^i, or see Chung and WilliamiP^ 
for a modern proof. We shall just remind the reader 
of this result, in its simplest terms. Suppose that W K 
is a Wiener measure over a rescaled Brownian motion 
\/2~K,W(t) and Wi/ 2 = W is the standard Wiener 
measure. Suppose that u(x, t) is any smooth vector field 
and W(t) is defined in terms of X(i) by integrating ffify: 



W(t) 



1 



/2k 



X(t)-x - / dsu(X{s),s) 
Jo 



Then Girsanov's theorem states that 

,-' I' 



Z>W K (X)exp 



M/ u(X,i)-dX- i / u 2 (K,t)dt 
' K \ Jo 2 J 

= X>W(W), (A.4) 



where 



/ u(X,t)-dX = lim Vu(X n _ 1 ,t n _ 1 )-(X n -X„_ 1 ) 

is the usual Ito stochastic integral 23 * 124 *. This is eas- 
ily seen to be equivalent to (A. 3 1 by expanding the 



square in the exponent of the latter and noting that 
DXexp f-j^j J2n=i l X « - x «-i| 2 ) converges to the 
Wiener measure 2?W K (X) in the continuum limit. 

Although the Euler-Maruyama scheme provides the 
simplest derivation of such path-integral formulas, other 
discretizations are possible and yield the same results. 
E.g., suppose that the trapezoidal rule is employed: 



X„ — X„ 



u(X„,£ ra ) + u(X„_i,t n _i j 



St 



^(W n - W n _r). (A.5) 



The sum in the exponent of (A. 3 1 is replaced by 

N 



i u (X n , t n ) + u(X n _i, t n _i) 
<5i ~ 2 



Expanding and taking the continuum limit, this con- 
verges formally to 



/o 
where 

l-T 



±\ 2 dt-2 / u(X,i)odX 



u 2 (X,t)dt, 



[ u(X,t)odX 
Jo 



N 



,. u(X n ,i n ) + u(X n _i,t n _i) 
hm > (X n -X n _i) 

iv— >-oo * — ' 7 



(A.6) 

is the Stratonovich stochastic integral. Since formally 

X>Q K (X) = PXexp f-i / Q T |X| 2 <ft) , it appears that 

we recover the same result as before but with the 
Stratonovich integral substituted for the Ito integral in 
equation (A.4 1. However, we have not yet computed the 



contribution of the Jacobian determinant to the change of 
variables. With the discretization ( A.5[ ) the Jacobian ma- 
trix is block lower-triangular, with diagonal blocks 

dX„ ,— ( 1 Su ; A -1 

^ = ^(i- 2 ^^(x„,uj 

= V2nexp (1st ^"(X,,, t n ) + 0(St 2 ] 
\2 ox 

Thus, using det (exp A) = exp(tr A), 

^exp^ T (V-u)(X,t)^ 

in the continuum limit N — > oo. The integral in the lat- 
ter exponent can be expressed in terms of the quadratic 



variation process 

(x,y) t = lim V(jf n -x„_i)(y n -y„_i). 

TV— voo * — * 

n=l 

Using d(Xi, Xj) t — 2n5ijdt, one obtains 

- | (V-u)(X, i) eft = i- £ d(u(X, t); X) . 

But the standard relation between Ito and Stratonovich 
integrals, 

/ u(X,i)odX-i [ d(u(X,i);X) = / u(X,t)-dX, 
Jo 2 j j 



then recovers (A.4 1 exactly as before 



For physicists, an illuminating derivation of the path- 
integral formulas can be based on Feynman's famous for- 
mula for the transition amplitude of a quantum, nonrel- 
ativistic, charged particle moving in a scalar potential V 
and in a magnetic field with vector potential A. Feyn- 
man's result was 



(x,i|x o ,0) 



c(t)= 



2?x exp 



ds C(s) 



<(0)=x \' l J0 

where C(t) is the classical Lagrangian 

C(t) = ^mlxl 2 + -A(x, f)-x - V(x, t) 
2 c 

and the amplitude satisfies the Schrodingcr equation 

ihd t ^ = — (-ihV - - A(x, t)\ 2 * + V(x, t)tf . 
2m V c / 
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See FeynmarPH, p. 376 or later treatment d 126 ' 127 1 . Taking 
imaginary time T = it, u 



rnc i 



and 



He 



V = -t^A 2 -i-^V-A, 
zmr 2mc 

yields the path-integral formula ([9| and converts the 
Schrodinger equation into the diffusion equation = 
kA\1/ — (u-V)^. This is very straightforward to check in 
the Coulomb gauge V-A = 0. For a general choice of 
gauge, note that the vector potential term in the classi- 
cal action C(s) yields a term in the exponent of Feyn- 
man's path-integral formula which must be interpreted 
as a Stratonovich integral: 

/ A(x, s)-x(s) ds = / A(x, s)o dx(s). 
Jo Jo 

This is implicit in Feynman's original derivation, who 
used t he m idpoint discretization to define the above 
integraP ^. Th is point has been carefully discussed 
elsewher e! 126 ! 127 !. If this Stratonovich integral is combined 
with the V-A term from the potential, one gets a net 
contribution to the action proportional to 



u(X,t)odX+- / d{u(X,t);X) = u(X,t)-dX, 
i) 2 J J 



where now t = is and 



/ u(X,t)-5x= lim Vu(X n ,t n )-(X n -X n _ 1 ) 
Jo N ^°°^i 

is the backward Ito integral. One obtains a result just like 
the Girsanov formula ( |A.4| but with the usual (forward) 
Ito integral replaced by a backward Ito integral. 

Feynman's result is correct. As we stressed in the text, 
our path-integral formulas Q and (111 correspond to 
solving the SDE ^ backward in time, e.g. with the 
backward Euler-Maruyama scheme 



X n - 
for L 



X n - u(X n ,t n )St + v^(W n _i - W„) (A.7) 



= t. 



5t. If we repeat the steps that led us to 



(A.3I, we now obtain 



V{Xi,...,X N ) oc 

N 



exp 



4k ^ 



n=l 



X„ 



St 



— - u(X„,t n ) 



(A. 



Integrating this density with respect to DX = J\ dX n 
and taking the continuous-time limit St — > 0, N — > oo, 
one formally obtains a Girsanov-type formula with the 
forward Ito integral replaced by a backward Ito integral. 
This is the rigorous version of our formulas ([9| and (11 1, 
which correspond exactly to Feynman's old result. 
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